//
// Created by tong on 9/27/24.
//

#include "GRDECLGrid.h"

#include <functional>
#include <iostream>
#include <sstream>
#include <fstream>
#include <algorithm>
#include <iomanip>

namespace{
    constexpr double eps = 2;
    inline bool sameEdge(const int a1, const int a2, const int b1, const int b2)
    {
        return (a1 == b1 && a2 == b2) || (a1 == b2 && a2 == b1);
    }

    inline bool coEdgeFace(const IntList &f1, const IntList &f2)
    {
        for (int i = 0; i < f1.size(); i++)
        {
            int a1 = f1[i];
            int a2 = i == (f1.size()-1) ? f1[0] : f1[i + 1];
            for (int j = 0; j < f2.size(); j++)
            {
                int b1 = f2[j];
                int b2 = j == (f2.size()-1) ? f2[0] : f2[j + 1];
                if (sameEdge(a1, a2, b1, b2))
                {
                    return true;
                }
            }
        }
        return false;
    }
    inline bool coPlaneFace(const Vec &f1n, const Vec &f2n)
    {
        Vec n1 = f1n / mag(f1n);
        Vec n2 = f2n / mag(f2n);
        return mag(n1 ^ n2) < 0.2;
    }
    inline bool coLineVec(const Vec &v1, const Vec &v2)
    {
        return VecAngle(v1, v2) < PI/18 || VecAngle(v1, v2) > PI*17/18;
    }

    template <typename T>
    T prevIdx(T idx, T size)
    {
        return idx == 0 ? size - 1 : idx - 1;
    }
    template <typename T>
    T nextIdx(T idx, T size)
    {
        return idx == size - 1 ? 0 : idx + 1;
    }
    double calcCellVolumeByNodes(
        const PointList &nodeCoord,
        const Point &cellCenter,
        const IntList &cellNodes)
    {
        double volume = 0.0;
        if (cellNodes.size() == 8)
        {
            volume += pyramidVolume(cellCenter,
                                    nodeCoord[cellNodes[0]],
                                    nodeCoord[cellNodes[1]],
                                    nodeCoord[cellNodes[2]],
                                    nodeCoord[cellNodes[3]]);
            volume += pyramidVolume(cellCenter,
                                    nodeCoord[cellNodes[4]],
                                    nodeCoord[cellNodes[5]],
                                    nodeCoord[cellNodes[6]],
                                    nodeCoord[cellNodes[7]]);
            volume += pyramidVolume(cellCenter,
                                    nodeCoord[cellNodes[0]],
                                    nodeCoord[cellNodes[1]],
                                    nodeCoord[cellNodes[5]],
                                    nodeCoord[cellNodes[4]]);
            volume += pyramidVolume(cellCenter,
                                    nodeCoord[cellNodes[1]],
                                    nodeCoord[cellNodes[2]],
                                    nodeCoord[cellNodes[6]],
                                    nodeCoord[cellNodes[5]]);
            volume += pyramidVolume(cellCenter,
                                    nodeCoord[cellNodes[2]],
                                    nodeCoord[cellNodes[3]],
                                    nodeCoord[cellNodes[7]],
                                    nodeCoord[cellNodes[6]]);   
            volume += pyramidVolume(cellCenter,
                                    nodeCoord[cellNodes[3]],
                                    nodeCoord[cellNodes[0]],
                                    nodeCoord[cellNodes[4]],
                                    nodeCoord[cellNodes[7]]);
        }
        else if(cellNodes.size() == 6)
        {
            volume += tetraVolume(cellCenter,
                                  nodeCoord[cellNodes[0]],
                                  nodeCoord[cellNodes[1]],
                                  nodeCoord[cellNodes[2]]);
            volume += tetraVolume(cellCenter,
                                  nodeCoord[cellNodes[0]],
                                  nodeCoord[cellNodes[2]],
                                  nodeCoord[cellNodes[3]]);
            volume += pyramidVolume(cellCenter,
                                    nodeCoord[cellNodes[0]],
                                    nodeCoord[cellNodes[1]],
                                    nodeCoord[cellNodes[4]],
                                    nodeCoord[cellNodes[3]]);
            volume += pyramidVolume(cellCenter,
                                    nodeCoord[cellNodes[1]],
                                    nodeCoord[cellNodes[2]],
                                    nodeCoord[cellNodes[5]],
                                    nodeCoord[cellNodes[4]]);
            volume += pyramidVolume(cellCenter,
                                    nodeCoord[cellNodes[2]],
                                    nodeCoord[cellNodes[0]],
                                    nodeCoord[cellNodes[3]],
                                    nodeCoord[cellNodes[5]]);
        }
        else if(cellNodes.size() == 4)
        {
            volume = tetraVolume(nodeCoord[cellNodes[0]],
                                 nodeCoord[cellNodes[1]],
                                 nodeCoord[cellNodes[2]],
                                 nodeCoord[cellNodes[3]]);
        }
        else if(cellNodes.size() ==5)
        {
            volume = pyramidVolume(nodeCoord[cellNodes[4]],
                                   nodeCoord[cellNodes[0]],
                                   nodeCoord[cellNodes[1]],
                                   nodeCoord[cellNodes[2]],
                                   nodeCoord[cellNodes[3]]);
        }
        if (volume < 0)
        {
            volume = 1e-13;
            std::cerr << "Negative volume!\n";
        }
        if (volume > 1e20)
        {
            volume = 0;
            std::cerr << "Too large volume!\n";
        }
        return volume;
    }
    
    void reorderCellNodes(const PointList& points,IntList& cellNodes){
        
        // std::sort(cellNodes.begin(),cellNodes.end(),[&points](int a,int b)
        // {
        //     // sort by z then y then x
        //     if ((points[a].z - points[b].z) > 1)
        //     {
        //         return true;
        //     }
        //     else if (std::abs(points[a].z - points[b].z) < 1){
        //         if ((points[a].y - points[b].y) > 1e-2)
        //         {
        //             return true;
        //         }
        //         else if (std::abs(points[a].y - points[b].y) < 1e-2)
        //         {
        //             if ((points[a].x - points[b].x) > 1e-2)
        //             {
        //                 return true;
        //             }
        //         }
        //     }
        //     return false;

        // });

    // sort by z
    std::sort(cellNodes.begin(), cellNodes.end(), [&points](int a, int b)
            { return points[a].z < points[b].z; });
    // sort by y and x
    std::sort(cellNodes.begin(), cellNodes.begin()+4, [&points](int a, int b)
            { if((points[a].y - points[b].y) < -1e-2)
            {return true;}
            else if(std::abs(points[a].y - points[b].y) < 1e-2){
                if((points[a].x - points[b].x) < 1e-2)
                {return true;}
            }
            return false; });
    std::sort(cellNodes.begin()+4, cellNodes.begin()+8, [&points](int a, int b)
              { if((points[a].y - points[b].y) < -1e-2)
            {return true;}
            else if(std::abs(points[a].y - points[b].y) < 1e-2){
                if((points[a].x - points[b].x) < -1e-2)
                {return true;}
            }
            return false; });        
    std::swap(cellNodes[2],cellNodes[3]);
    std::swap(cellNodes[6],cellNodes[7]);
    }
    
}
double GRDECLGrid::tol_ = 2.2204e-6;
bool GRDECLGrid::perfFlag_ = false;
bool GRDECLGrid::pinchActive_ = true;
bool GRDECLGrid::writeDebugFiles_ = false;
GRDECLGrid::CellCenterType GRDECLGrid::cellCenterType_ = CellCenterType::Barycenter;
double GRDECLGrid::extendZMax_ = 1200;
double GRDECLGrid::extendZMin_ = -1000;
double GRDECLGrid::extendXFactor_ = 0.1;
double GRDECLGrid::extendYFactor_ = 0.1;
int GRDECLGrid::XExtendNum_ = 5;
int GRDECLGrid::YExtendNum_ = 5;
int GRDECLGrid::ZExtendNum_ = 4;

inline Point interpolatePillar(const double *coord, const double zCoord)
{
    const double dz = coord[5] - coord[2];
    double a;
    Point pt(0, 0, zCoord);
    if (std::fabs(dz > 0))
    {
        a = (pt.z - coord[2]) / dz;
    }
    else
    {
        a = 0;
    }
    pt.x = (1.0 - a) * coord[0] + a * coord[3];
    pt.y = (1.0 - a) * coord[1] + a * coord[4];
    return pt;
}

inline bool lineIntersect(const int a1, const int a2, const int b1, const int b2)
{
    return (((a1 > b1) && (a2 < b2)) || ((a1 < b1) && (a2 > b2)));
}

inline bool faceIntersect(
    const int a10, const int a20, const int a11, const int a21,
    const int b10, const int b20, const int b11, const int b21)
{

    return max(a10, b10) < min(a11, b11) ||
           max(a20, b20) < min(a21, b21) ||
           lineIntersect(a10, a20, b10, b20) ||
           lineIntersect(a11, a21, b11, b21);
}

inline bool validCellFace(const int a10, const int a11, const int b10, const int b11)
{
    return (!((a10 == INT32_MIN) && (b10 == INT32_MIN)) && !((a11 == INT32_MAX) && (b11 == INT32_MAX)));
}

inline GRDECLGrid::Face computeFaceTopology(
    const int a10, const int a20, const int a11, const int a21,
    const int b10, const int b20, const int b11, const int b21, IntersectIndex intersect)
{
    array<int, 8> mask = {-1};
    int k;
    GRDECLGrid::Face f;
    /* Which pillar points should we use */
    if (a11 > b11)
    {
        mask[0] = b11;
    }
    else
    {
        mask[0] = a11;
    }
    if (a21 > b21)
    {
        mask[2] = b21;
    }
    else
    {
        mask[2] = a21;
    }
    if (a20 > b20)
    {
        mask[4] = a20;
    }
    else
    {
        mask[4] = b20;
    }
    if (a10 > b10)
    {
        mask[6] = a10;
    }
    else
    {
        mask[6] = b10;
    }
    /* Illegal situations */

    if (mask[0] == mask[2] ||
        mask[0] == mask[4] ||
        mask[2] == mask[6] ||
        mask[4] == mask[6])
    {
        std::cerr << "Illegal Partial pinch!\n";
        abort();
    }

    /* Partial pinch of face */
    if (mask[0] == mask[6])
    {
        mask[6] = -1;
    }
    if (mask[2] == mask[4])
    {
        mask[4] = -1;
    }
    /* Get shape of face: */
    /*   each new intersection will be part of the new face, */
    /*   but not all pillar points. This is encoded in mask. */
    mask[1] = intersect[3]; /* top-top */
    mask[3] = -1;
    mask[5] = intersect[0]; /* bottom-bottom*/
    mask[7] = -1;
    /* bottom-top */
    if (intersect[1] != -1)
    {
        if (a10 > b11)
        { /* intersection[1] left of (any) intersection[0] */
            mask[0] = -1;
            mask[6] = -1;
            mask[7] = intersect[1];
        }
        else
        {
            mask[2] = -1;
            mask[4] = -1;
            mask[3] = intersect[1];
        }
    }
    /* top-bottom */
    if (intersect[2] != -1)
    {
        if (a11 < b10)
        { /* intersection[2] left of (any) intersection[3] */
            mask[0] = -1;
            mask[6] = -1;
            mask[7] = intersect[2];
        }
        else
        {
            mask[2] = -1;
            mask[4] = -1;
            mask[3] = intersect[2];
        }
    }
    for (k = 7; k >= 0; --k)
    {
        if (mask[k] != -1)
        {
            f.push_back(mask[k]);
        }
    }
    return f;
}

GRDECLGrid::GRDECLGrid(const GRDECLRawData &rawGRDECL)
    : rawGRDECL_(rawGRDECL),
      nRawCell_(rawGRDECL.dimI * rawGRDECL.dimJ * rawGRDECL.dimK),
      zcornSign_(getZCORNSign()),
      CoordSysType_(getCoordSysType()),
      nCell_(0),
      globalCellLocalIndex_(nRawCell_, -1),
      cornerMap_(rawGRDECL.dimJ + 1, CornerListList(rawGRDECL.dimI + 1)),
      posPointIndex_(rawGRDECL.dimJ * 2, IntListList(rawGRDECL.dimI * 2))
{
    getNodeCoord();
    getVerticalFace();
    getHorizontalFaces();
    calcIntersectPoints();
    reorderGRDECLData();

    calcPolyGridData();

    if (rawGRDECL_.YoungsMod_.size() == size_t(nRawCell_) &&
        rawGRDECL_.Poisson_.size() == size_t(nRawCell_) &&
        rawGRDECL_.RockDensity_.size() == size_t(nRawCell_) &&
        rawGRDECL_.Porosity_.size() == size_t(nRawCell_))
    {
        YoungsMod_.resize(nCell_);
        Poisson_.resize(nCell_);
        RockDens_.resize(nCell_);
        Poro_.resize(nCell_);
        // for younmod and poisson
        for (int i = 0; i < nCell_; i++)
        {
            YoungsMod_[i] = rawGRDECL_.YoungsMod_[localCellGlobalIndex_[i]];
            Poisson_[i] = rawGRDECL_.Poisson_[localCellGlobalIndex_[i]];
            Poro_[i] = rawGRDECL_.Porosity_[localCellGlobalIndex_[i]];
            RockDens_[i] = rawGRDECL_.RockDensity_[localCellGlobalIndex_[i]];
        }
        std::cout << "get YoungsMod_ and Poisson_\n";
    }
    else if (rawGRDECL_.YoungsMod_.empty() && rawGRDECL_.Poisson_.empty())
    {
        std::cout << "YoungsMod_ and Poisson_ is empty\n";
    }
    else
    {
        std::cerr << "YoungsMod_ and Poisson_ size is not match\n";
        std::abort();
    }
}

int GRDECLGrid::getZCORNSign() const
{
    /* Ensure that zcorn (depth) is strictly nondecreasing in
       the k-direction.
       1) if  z(i,j,k) <= z(i,j,k+1) for all (i,j,k), return 1.0
       2) if -z(i,j,k) <=-z(i,j,k+1) for all (i,j,k), return -1.0
       3) if (1) and (2) fails, return -1.0, and set *error = 1.
    */
    clock_t startT = std::clock();
    const DoubleList &zcorn = rawGRDECL_.zcoord;
    const BoolList &actnum = rawGRDECL_.actnum;
    const int nx = rawGRDECL_.dimI;
    const int ny = rawGRDECL_.dimJ;
    const int nz = rawGRDECL_.dimK;
    // std::cout<<rawGRDECL_.actnum.size()<<std::endl;
    assert(!actnum.empty());
    int sign;
    bool nondecreas = true;
    for (sign = 1; sign > -2; sign = sign - 2)
    {
        for (int j = 0; j < 2 * ny; ++j)
        {
            for (int i = 0; i < 2 * nx; ++i)
            {
                for (int k = 0; k < 2 * nz - 1; ++k)
                {
                    double z1 = sign * zcorn[i + 2 * nx * (j + 2 * ny * (k))];
                    double z2 = sign * zcorn[i + 2 * nx * (j + 2 * ny * (k + 1))];
                    int c1 = i / 2 + nx * (j / 2 + ny * (k / 2));
                    int c2 = i / 2 + nx * (j / 2 + ny * ((k + 1) / 2));
                    assert(c1 < (nx * ny * nz));
                    assert(c2 < (nx * ny * nz));
                    if ((actnum[c1] && actnum[c2]) && (z2 < z1))
                    {
                        std::cout << "\nZCORN should be strictly nondecreasing along pillars! "
                                     "Attempt to reverse sign in ZCORN\n";
                        if (sign == -1)
                            nondecreas = false;
                        goto end;
                    }
                }
            }
        }
        if (sign == 1 && nondecreas == true)
        {
            break;
        }
    end:
        break;
    }
    if (nondecreas == false)
    {
        std::cerr << "Failed to reverse sign in ZCORN\n";
        std::abort();
    }

    clock_t endT1 = std::clock();
    if (perfFlag_)
    {
        std::cout << "getZCORNSign() Time: " << double(endT1 - startT) / CLOCKS_PER_SEC << " s" << std::endl;
    }
    return sign;
}

GRDECLGrid::CoordinateSystemType GRDECLGrid::getCoordSysType() const
{
    const size_t nx = rawGRDECL_.dimI;
    const size_t ny = rawGRDECL_.dimJ;
    const size_t nc = nRawCell_;
    const BoolList &actnum = rawGRDECL_.actnum;
    size_t c = 0;
    size_t off[8];

    off[0] = 0;
    off[1] = off[0] + 1;
    off[2] = off[0] + (2 * nx);
    off[3] = off[2] + 1;
    off[4] = off[0] + ((2 * nx) * (2 * ny));
    off[5] = off[4] + 1;
    off[6] = off[4] + (2 * nx);
    off[7] = off[6] + 1;
    double dz = 0.0;
    const double *pt_coord = rawGRDECL_.coord.data();
    size_t origin = 0;
    size_t imax = (nx + 0) * 1 * (2 * 3);
    size_t jmax = (nx + 1) * (ny + 0) * (2 * 3);
    double dx[2], dy[2];
    dx[0] = pt_coord[imax + 0] - pt_coord[origin + 0];
    dy[0] = pt_coord[imax + 1] - pt_coord[origin + 1];
    dx[1] = pt_coord[jmax + 0] - pt_coord[origin + 0];
    dy[1] = pt_coord[jmax + 1] - pt_coord[origin + 1];

    bool active, searching;
    do
    {
        active = (actnum[c] != 0);
        if (active)
        {
            dz = vertSize(c, off);
        }
        searching = !(active && (fabs(dz) > 0.0));
        c += 1;
    } while (searching && (c < nc));

    assert(!searching); /* active && (fabs(dz) > 0) */

    /* Compute vector triple product to distinguish left-handed (<0)
     * from right-handed (>0) coordinate systems. */
    double triple = zcornSign_ * dz * (dx[0] * dy[1] - dx[1] * dy[0]);

    if (triple > 0.0)
    {
        return CoordinateSystemType::RightHanded;
    }
    else if (triple < 0.0)
    {
        return CoordinateSystemType::LeftHanded;
    }
    else
    {
        std::cerr << "\nCoordinateSystemType::Inconclusive!\n";
        return CoordinateSystemType::Inconclusive;
    }
}

void GRDECLGrid::getNodeCoord()
{
    clock_t startT = std::clock();
    const int nZPerLayer = 4 * rawGRDECL_.dimI * rawGRDECL_.dimJ;
    for (int zidx = 0; size_t(zidx) < rawGRDECL_.zcoord.size(); zidx++)
    {
        const int nk = zidx / nZPerLayer;
        const int K = nk / 2;
        const int posK = (nk + 1) / 2;
        const int layerZidx = zidx % nZPerLayer;
        const int nJ = layerZidx / (2 * rawGRDECL_.dimI);
        const int J = nJ / 2;
        const int pillarJ = (nJ + 1) / 2;
        const int nI = layerZidx % (2 * rawGRDECL_.dimI);
        const int I = nI / 2;
        const int pillarI = (nI + 1) / 2;
        const int cellIndex = cellLinearIndex(I, J, K);
        cornerMap_[pillarJ][pillarI].push_back(CornerPoint(rawGRDECL_.zcoord[zidx], cellIndex, nI, nJ, posK, rawGRDECL_.actnum[cellIndex]));
    }

    clock_t endT1 = std::clock();
    if (perfFlag_)
    {
        std::cout << "recon getNodeCoord insert point:: " << double(endT1 - startT) / CLOCKS_PER_SEC << " s" << std::endl;
    }

    // int pNum = 0;
    for (int pillarJ = 0; pillarJ < int(cornerMap_.size()); pillarJ++)
    {
        for (int pillarI = 0; pillarI < int(cornerMap_[pillarJ].size()); pillarI++)
        {
            // int pillarIndex = pillarJ + pillarI*(rawGRDECL_.dimI+1);
            const CornerList &cpList = cornerMap_[pillarJ][pillarI];
            // std::cout<<"pillarI:" <<pillarI<<" pillarJ:"<<pillarJ<<" size():"<<cpList.size()<<"\n";
            // if(cpList.empty())
            //     continue;
            const double *pCoord = &rawGRDECL_.coord[6 * ((rawGRDECL_.dimI + 1) * pillarJ + pillarI)];
            DoubleList pillarZcoord;
            // add deepth coord
            for (int nz = 0; nz < int(cpList.size()); nz++)
            {
                const CornerPoint &cp = cpList[nz];
                if (cp.active)
                {
                    pillarZcoord.push_back(cp.zcoord);
                }
            }
            if (pillarZcoord.empty())
            {
                for (int nz = 0; nz < int(cpList.size()); nz++)
                {   
                    // get all points
                    const CornerPoint &cp = cpList[nz];
                    if (posPointIndex_[cp.posJ][cp.posI].size() == 0)
                    {
                        posPointIndex_[cp.posJ][cp.posI].emplace_back(INT32_MIN);
                    }
                    posPointIndex_[cp.posJ][cp.posI].emplace_back(posPointIndex_[cp.posJ][cp.posI].back());
                    if (posPointIndex_[cp.posJ][cp.posI].size() == size_t(2 * rawGRDECL_.dimK + 1))
                    {
                        posPointIndex_[cp.posJ][cp.posI].emplace_back(INT32_MAX);
                    }
                }
                continue;
            }
            IntList orderIndex = sortIndex(pillarZcoord);
            IntList activeNodeIndex;
            double val = pillarZcoord[orderIndex[0]]; /* Keep first value */
            // const CornerPoint& cp1=cpList[orderIndex[0]];
            nodeCoord_.emplace_back(interpolatePillar(pCoord, val));
            activeNodeIndex.push_back(int(nodeCoord_.size()) - 1);
            // posPointIndex_[cp1.posJ][cp1.posI].emplace_back(nodeCoord_.size() - 1);
            for (int zci = 1; zci < int(orderIndex.size()); zci++)
            {
                if (val + tol_ < pillarZcoord[orderIndex[zci]])
                {
                    val = pillarZcoord[orderIndex[zci]];
                    const Point pt(interpolatePillar(pCoord, val));
                    nodeCoord_.emplace_back(pt);
                }
                activeNodeIndex.push_back(int(nodeCoord_.size()) - 1);
            }
            int actIdx = 0;
            for (int nz = 0; nz < int(cpList.size()); nz++)
            {
                const CornerPoint &cp = cpList[nz];
                if (posPointIndex_[cp.posJ][cp.posI].size() == 0)
                {
                    posPointIndex_[cp.posJ][cp.posI].push_back(INT32_MIN);
                }
                if (cp.active)
                {
                    posPointIndex_[cp.posJ][cp.posI].push_back(activeNodeIndex[actIdx++]);
                }
                else
                {
                    posPointIndex_[cp.posJ][cp.posI].push_back(posPointIndex_[cp.posJ][cp.posI].back());
                }
                if (posPointIndex_[cp.posJ][cp.posI].size() == size_t(2 * rawGRDECL_.dimK + 1))
                {
                    posPointIndex_[cp.posJ][cp.posI].push_back(INT32_MAX);
                }
            }
        }
    }
    clock_t endT2 = std::clock();
    if (perfFlag_)
    {
        std::cout << "recon getNodeCoord calc unique point: " << double(endT2 - endT1) / CLOCKS_PER_SEC << " s" << std::endl;
    }

    if (CoordSysType_ == CoordinateSystemType::LeftHanded)
    {
        for (size_t i = 0; i < nodeCoord_.size(); i++)
        {
            nodeCoord_[i].y = -nodeCoord_[i].y;
        }
    }

    if (writeDebugFiles_)
    {
        std::ofstream outfile("reconNodes.txt", std::ios::out | std::ios::trunc);
        outfile << nodeCoord_.size() << "\n";
        for (int i = 0; i < int(nodeCoord_.size()); i++)
        {
            outfile << nodeCoord_[i].x << ", " << nodeCoord_[i].y << ", " << nodeCoord_[i].z << "\n";
        }
        outfile << "\\ \n";

        std::ofstream outfile2("reconPos.txt", std::ios::out | std::ios::trunc);
        int numPl = 0;
        for (int posJ = 0; posJ < int(posPointIndex_.size()); posJ++)
        {
            for (int posI = 0; posI < int(posPointIndex_[posJ].size()); posI++)
            {
                for (int i = 0; i < int(posPointIndex_[posJ][posI].size()); i++)
                {
                    outfile2 << posPointIndex_[posJ][posI][i] << "\n";
                    numPl++;
                }
            }
        }
        outfile2 << "\\ " << numPl << "\n";
    }
}

void GRDECLGrid::getVerticalFace()
{

    clock_t startT = std::clock();
    const int dimI = rawGRDECL_.dimI;
    const int dimJ = rawGRDECL_.dimJ;
    // const int dimK = rawGRDECL_.dimK;
    //  for internal face in I direction
    // const int IFaceStart = faces_.size();
    for (int j = 0; j < dimJ; j++)
    { // 0-dimJ
        for (int i = 0; i < dimI + 1; i++)
        { // 0-dimI+1
            const IntList &a1 = posPointIndex_[2 * j][max(0, 2 * i - 1)];
            const IntList &a2 = posPointIndex_[2 * j + 1][max(0, 2 * i - 1)];
            const IntList &b1 = posPointIndex_[2 * j][min(2 * i, dimI * 2 - 1)];
            const IntList &b2 = posPointIndex_[2 * j + 1][min(2 * i, dimI * 2 - 1)];
            IntList czIndex = findVerticalFace(a1, a2, b1, b2);
            assert((czIndex.size() % 2) == 0);
            for (int dk = 0; dk < int(czIndex.size()); dk += 2)
            {
                if (i == 0)
                {
                    czIndex[dk] = -1;
                    IMinFace_.insert(int(faceTags_.size()));
                }
                if (i == dimI)
                {
                    czIndex[dk + 1] = -1;
                    IMaxFace_.insert(int(faceTags_.size()));
                }
                faceNeighbors_.push_back(czIndex[dk] == -1 ? -1 : cellLinearIndex(i - 1, j, czIndex[dk]));
                faceNeighbors_.push_back(czIndex[dk + 1] == -1 ? -1 : cellLinearIndex(i, j, czIndex[dk + 1]));
                faceTags_.push_back(FaceType::IFace);
            }
        }
    }
    // faceTags_.resize(faceNeighbors_.size(),FaceType::InternalI);
    //  for internal face in J direction
    for (int j = 0; j < dimJ + 1; j++)
    { // 0-J+1
        for (int i = 0; i < dimI; i++)
        { // 0-I
            const IntList &a1 = posPointIndex_[max(0, 2 * j - 1)][2 * i + 1];
            const IntList &a2 = posPointIndex_[max(0, 2 * j - 1)][2 * i];
            const IntList &b1 = posPointIndex_[min(2 * j, dimJ * 2 - 1)][2 * i + 1];
            const IntList &b2 = posPointIndex_[min(2 * j, dimJ * 2 - 1)][2 * i];
            IntList czIndex = findVerticalFace(a1, a2, b1, b2);
            assert((czIndex.size() % 2) == 0);
            for (int dk = 0; dk < int(czIndex.size()); dk += 2)
            {
                if (j == 0)
                {
                    czIndex[dk] = -1;
                    JMinFace_.insert(int(faceTags_.size()));
                }
                if (j == dimJ)
                {
                    czIndex[dk + 1] = -1;
                    JMaxFace_.insert(int(faceTags_.size()));
                }
                faceNeighbors_.push_back(czIndex[dk] == -1 ? -1 : cellLinearIndex(i, j - 1, czIndex[dk]));
                faceNeighbors_.push_back(czIndex[dk + 1] == -1 ? -1 : cellLinearIndex(i, j, czIndex[dk + 1]));
                faceTags_.push_back(FaceType::JFace);
            }
        }
    }
    clock_t endT1 = std::clock();
    if (perfFlag_)
    {
        std::cout << "recon getVerticalFace: " << double(endT1 - startT) / CLOCKS_PER_SEC << " s" << std::endl;
    }

    if (writeDebugFiles_)
    {
        std::ofstream outfile("reconVerticalFaces.txt", std::ios::out | std::ios::trunc);
        outfile << faces_.size() << "\n";
        for (int i = 0; i < int(faces_.size()); i++)
        {
            for (int pi = 0; pi < int(faces_[i].size()); pi++)
                outfile << faces_[i][pi] << " ";
            outfile << "\n";
        }
    }
    if (faceTags_.size() != (faces_.size()))
    {
        std::cerr << "faceTags_ size is not equal to faceNeighbors_ size\n";
        std::abort();
    }
}

IntList GRDECLGrid::findVerticalFace(
    const IntList &a1, 
    const IntList &a2, 
    const IntList &b1, 
    const IntList &b2)
{
    IntList cZIndex;
    const int n = 2 * rawGRDECL_.dimK + 2;
    IntList interTop(n,-1);
    IntList interBottom(n,-1);
    IntersectIndex intersect;
    std::fill(intersect.begin(), intersect.end(), -1);
    int mj1 = 0;
    int mj2 = 0;
    int j = 0;
    for (int i = 0; i < n - 1; ++i)
    {
        /* pinched a-cell */
        if ((a1[i] == a1[i + 1]) &&
            (a2[i] == a2[i + 1]))
        {
            continue;
        }
        while ((j < n - 1) &&
               ((b1[j] < a1[i + 1]) ||
                (b2[j] < a2[i + 1])))
        {
            /* pinched b-cell */
            if ((b1[j] == b1[j + 1]) &&
                (b2[j] == b2[j + 1]))
            {
                interTop[j + 1] = interTop[j];
                j++;
                continue;
            }
            // face a(i,i+1) and face b(j,j+1) have nonzero intersection
            if (faceIntersect(a1[i], a2[i], a1[i + 1], a2[i + 1],
                              b1[j], b2[j], b1[j + 1], b2[j + 1]))
            {
                /* Completely matching faces */
                if (((a1[i] == b1[j]) && (a1[i + 1] == b1[j + 1])) &&
                    ((a2[i] == b2[j]) && (a2[i + 1] == b2[j + 1])))
                {
                    if (validCellFace(a1[i], a1[i + 1], b1[j], b1[j + 1]))
                    {
                        int cell_a = i % 2 != 0 ? (i - 1) / 2 : -1;
                        int cell_b = j % 2 != 0 ? (j - 1) / 2 : -1;
                        if (cell_a != -1 || cell_b != -1)
                        {
                            cZIndex.push_back(cell_a);
                            cZIndex.push_back(cell_b);
                            Face tempFace;
                            tempFace.push_back(a1[i]);
                            tempFace.push_back(a2[i]);
                            /* avoid duplicating nodes in pinched faces  */
                            if (a2[i + 1] != a2[i])
                            {
                                tempFace.push_back(a2[i + 1]);
                            }
                            if (a1[i + 1] != a1[i])
                            {
                                tempFace.push_back(a1[i + 1]);
                            }
                            // out->face_ptr[++out->number_of_faces] = f - out->face_nodes;
                            faces_.emplace_back(tempFace);
                        } // end for  cell_a != -1 || cell_b != -1
                    } // end for validCellFace(a1[i], a1[i + 1], b1[j], b1[j + 1])
                } // end for Completely matching faces
                /* Non-matching faces */
                else
                {
                    /* Find new intersection */
                    if (lineIntersect(a1[i + 1], a2[i + 1], b1[j + 1], b2[j + 1]))
                    {
                        /* store point numbers of intersecting lines */
                        interTop[j + 1] = int(nodeCoord_.size()) + int(intersectIndices_.size());
                        IntersectIndex newInter;
                        newInter[0] = a1[i + 1];
                        newInter[1] = a2[i + 1];
                        newInter[2] = b1[j + 1];
                        newInter[3] = b2[j + 1];
                        intersectIndices_.push_back(newInter);
                    }
                    else
                    {
                        interTop[j + 1] = -1;
                    }
                    /* Update intersection record */
                    intersect[0] = interBottom[j];     /* i   x j   */
                    intersect[1] = interBottom[j + 1]; /* i   x j+1 */
                    intersect[2] = interTop[j];        /* i+1 x j   */
                    intersect[3] = interTop[j + 1];    /* i+1 x j+1 */
                    /* Add face to list of faces if no INT_MIN or
                     * INT_MAX appear in a or b. */
                    if (validCellFace(a1[i], a1[i + 1], b1[j], b1[j + 1]))
                    {
                        /*Even indices refer to space between cells,
                        odd indices refer to cells*/
                        int cell_a = i % 2 != 0 ? (i - 1) / 2 : -1;
                        int cell_b = j % 2 != 0 ? (j - 1) / 2 : -1;
                        if (cell_a != -1 || cell_b != -1)
                        {
                            cZIndex.push_back(cell_a);
                            cZIndex.push_back(cell_b);
                            // computeFaceTopology(a1+i, a2+i, b1+j, b2+j, intersect, f);
                            if(&a1[i]!=&b1[j] && &a2[i]!=&b2[j]){
                                faultFaces_.push_back(int(faces_.size()));
                            }
                            faces_.emplace_back(computeFaceTopology(a1[i], a2[i], a1[i + 1], a2[i + 1],
                                                                    b1[j], b2[j], b1[j + 1], b2[j + 1], intersect));
                            
                        } // end for cell_a != -1 || cell_b != -1
                    } // end for validCellFace(a1[i], a1[i + 1], b1[j], b1[j + 1]
                } // end for Non-matching faces
            } // end for faceIntersect
            if (b1[j] < a1[i + 1])
            {
                mj1 = j;
            }
            if (b2[j] < a2[i + 1])
            {
                mj2 = j;
            }
            j++;
        } // end for while face bj
        std::copy(interBottom.begin(), interBottom.end(), interTop.begin());
        //interBottom = interTop;
        /* Zero out the "new" itop */
        std::fill(interTop.begin(), interTop.end(), -1);
        //interTop = {-1};
        j = min(mj1, mj2);
    } // end for face ai
    return cZIndex;
}

void GRDECLGrid::getHorizontalFaces()
{
    clock_t startT = std::clock();
    const int dimI = rawGRDECL_.dimI;
    const int dimJ = rawGRDECL_.dimJ;
    const int dimK = rawGRDECL_.dimK;

    int prevCell = -1;
    int thisCell = -1;

    int cellCounter = 0;
    for (int j = 0; j < dimJ; j++)
    {
        for (int i = 0; i < dimI; i++)
        {
            const IntList &c0 = posPointIndex_[2 * j][2 * i];
            const IntList &c1 = posPointIndex_[2 * j + 1][2 * i];
            const IntList &c2 = posPointIndex_[2 * j][2 * i + 1];
            const IntList &c3 = posPointIndex_[2 * j + 1][2 * i + 1];

            prevCell = -1;

            for (int k = 1; k < 2 * dimK + 1; k++)
            {
                const int globalCellIdx = cellLinearIndex(i, j, (k - 1) / 2);

                if (c0[k] == c0[k + 1] && c1[k] == c1[k + 1] &&
                    c2[k] == c2[k + 1] && c3[k] == c3[k + 1])
                {
                    if (k % 2)
                    {
                        globalCellLocalIndex_[globalCellIdx] = -1;
                    }
                }
                else
                {
                    if (k % 2)
                    {
                        thisCell = globalCellIdx;
                        if (!pinchActive_ &&
                            !isVerticalNeighbor(thisCell, prevCell) &&
                            prevCell != -1)
                        {
                            assert(faces_.size() > 0);
                            if (faceNeighbors_.back() == -1)
                            {
                                assert(faceNeighbors_[faceNeighbors_.size() - 2] == faceNeighbors_.back());
                            }
                            else
                            {
                                assert(faceNeighbors_.back() == prevCell);
                                if (k == 1)
                                {
                                    KMinFace_.insert(int(faceTags_.size()));
                                }
                                else if (k == 2 * dimK)
                                {
                                    KMaxFace_.insert(int(faceTags_.size()));
                                }
                                faceKValue_.emplace(int(faceTags_.size()), k / 2);
                                faces_.emplace_back(Face{c0[k], c2[k], c3[k], c1[k]});
                                faceTags_.push_back(FaceType::KFace);
                                faceNeighbors_.push_back(prevCell);
                                faceNeighbors_.push_back(-1);
                            }
                        }
                        if (k == 1)
                        {
                            KMinFace_.insert(int(faceTags_.size()));
                        }
                        else if (k == 2 * dimK)
                        {
                            KMaxFace_.insert(int(faceTags_.size()));
                        }
                        faceKValue_.emplace(int(faceTags_.size()), k / 2);
                        faces_.emplace_back(Face{c0[k], c2[k], c3[k], c1[k]});
                        faceTags_.push_back(FaceType::KFace);
                        faceNeighbors_.push_back(
                            (pinchActive_ || isVerticalNeighbor(thisCell, prevCell)) ? prevCell : -1);
                        faceNeighbors_.push_back(thisCell);
                        prevCell = thisCell;
                        globalCellLocalIndex_[thisCell] = cellCounter;
                        ++cellCounter;
                    }
                    else
                    {
                        if (prevCell != -1)
                        {   
                            if (k == 1)
                            {
                                KMinFace_.insert(int(faceTags_.size()));
                            }
                            else if (k == 2 * dimK)
                            {
                                KMaxFace_.insert(int(faceTags_.size()));
                            }
                            faceKValue_.emplace(int(faceTags_.size()), k / 2);
                            faces_.emplace_back(Face{c0[k], c2[k], c3[k], c1[k]});
                            faceTags_.push_back(FaceType::KFace);
                            faceNeighbors_.push_back(prevCell);
                            faceNeighbors_.push_back(-1);
                            prevCell = -1;
                        }
                    }

                }
            }
        }
    }
    nCell_ = cellCounter;

    if (faceTags_.size() != (faces_.size()))
    {
        std::cerr << "faceTags_ size is not equal to faceNeighbors_ size\n";
        std::abort();
    }

    clock_t endT1 = std::clock();
    if (perfFlag_)
    {
        std::cout << "recon getHorizontalFaces: " << double(endT1 - startT) / CLOCKS_PER_SEC << " s" << std::endl;
    }

    if (writeDebugFiles_)
    {
        std::ofstream outfile("reconAllFaces.txt", std::ios::out | std::ios::trunc);
        outfile << faces_.size() << "\n";
        for (int i = 0; i < int(faces_.size()); i++)
        {
            for (int pi = 0; pi < int(faces_[i].size()); pi++)
                outfile << faces_[i][pi] << " ";
            outfile << "\n";
        }
    }
}

void GRDECLGrid::calcIntersectPoints()
{
    clock_t startT = std::clock();
    for (size_t i = 0; i < intersectIndices_.size(); i++)
    {
        const IntersectIndex &itsc = intersectIndices_[i];
        const Point pt0 = nodeCoord_[itsc[0]];
        const Point pt1 = nodeCoord_[itsc[1]];
        const Point pt2 = nodeCoord_[itsc[2]];
        const Point pt3 = nodeCoord_[itsc[3]];
        double a;
        if (fabs((pt1.z - pt0.z) - (pt3.z - pt2.z)) > 0.0)
        {
            a = (pt2.z - pt0.z) / ((pt1.z - pt0.z) - (pt3.z - pt2.z));
        }
        else
        {
            a = 0;
        }
        const double z = pt0.z * (1.0 - a) + pt1.z * a;
        double b1 = (pt2.z - z) / (pt2.z - pt0.z);
        double b2 = (z - pt0.z) / (pt2.z - pt0.z);
        const double x1 = pt0.x * b1 + pt2.x * b2;
        const double y1 = pt0.y * b1 + pt2.y * b2;
        b1 = (z - pt3.z) / (pt1.z - pt3.z);
        b2 = (pt1.z - z) / (pt1.z - pt3.z);
        const double x2 = pt1.x * b1 + pt3.x * b2;
        const double y2 = pt1.y * b1 + pt3.y * b2;
        nodeCoord_.emplace_back(Point(x1 * (1.0 - a) + x2 * a, y1 * (1.0 - a) + y2 * a, z));
    }
    clock_t endT1 = std::clock();
    if (perfFlag_)
    {
        std::cout << "recon calcIntersectPoints: " << double(endT1 - startT) / CLOCKS_PER_SEC << " s" << std::endl;
    }

    if (writeDebugFiles_)
    {
        std::ofstream outfile("allReconNodes.txt", std::ios::out | std::ios::trunc);
        outfile << nodeCoord_.size() << "\n";
        for (int i = 0; i < int(nodeCoord_.size()); i++)
        {
            outfile << nodeCoord_[i].x << ", " << nodeCoord_[i].y << ", " << nodeCoord_[i].z << "\n";
        }
        outfile << "\\ \n";
    }
}

void GRDECLGrid::reorderGRDECLData()
{
    // shirnk to fit
    nodeCoord_.shrink_to_fit();
    faceNeighbors_.shrink_to_fit();
    faceTags_.shrink_to_fit();
    faces_.shrink_to_fit();
    for (size_t i = 0; i < faces_.size(); i++)
    {
        faces_[i].shrink_to_fit();
    }

    clock_t startT = std::clock();
    localCellGlobalIndex_.resize(nRawCell_, -1);
    int cellNum = 0;
    for (size_t i = 0; i < globalCellLocalIndex_.size(); i++)
    {
        if (globalCellLocalIndex_[i] != -1)
        {
            localCellGlobalIndex_[cellNum] = int(i);
            globalCellLocalIndex_[i] = cellNum;
            ++cellNum;
        }
    }
    localCellGlobalIndex_.resize(size_t(cellNum));

    for (size_t i = 0; i < faceNeighbors_.size(); i++)
    {
        if (faceNeighbors_[i] != -1)
        {
            faceNeighbors_[i] = globalCellLocalIndex_[faceNeighbors_[i]];
        }
    }

    if (CoordSysType_ == CoordinateSystemType::LeftHanded)
    {
        for (size_t i = 0; i < nodeCoord_.size(); i++)
        {
            nodeCoord_[i].y = -nodeCoord_[i].y;
        }
    }

    if (zcornSign_ == -1)
    {
        for (size_t i = 0; i < nodeCoord_.size(); i++)
        {
            nodeCoord_[i].z = -nodeCoord_[i].z;
        }
    }

    if ((CoordSysType_ == CoordinateSystemType::LeftHanded) ^ (zcornSign_ == -1))
    {
        for (size_t i = 0; i < faces_.size(); i++)
        {
            std::reverse(faces_[i].begin(), faces_[i].end());
        }
    }

    clock_t endT1 = std::clock();
    if (perfFlag_)
    {
        std::cout << "reorderGRDECLData Time: " << double(endT1 - startT) / CLOCKS_PER_SEC << " s" << std::endl;
    }

    if (writeDebugFiles_)
    {
        std::ofstream outfile("reconNeighbor.txt", std::ios::out | std::ios::trunc);
        outfile << faceNeighbors_.size() << "\n";
        for (size_t i = 0; i < faceNeighbors_.size(); i += 2)
        {
            outfile << faceNeighbors_[i] << " " << faceNeighbors_[i + 1] << "\n";
        }
        outfile << "\\ \n";

        std::ofstream outfile1("reconCellGlobalIndex.txt", std::ios::out | std::ios::trunc);
        outfile1 << localCellGlobalIndex_.size() << "\n";
        for (size_t i = 0; i < localCellGlobalIndex_.size(); i++)
        {
            outfile1 << localCellGlobalIndex_[i] << "\n";
        }
        outfile1 << "\\ \n";
    }
}

void GRDECLGrid::calcPolyGridData()
{

    clock_t startT = std::clock();

    // calc faceCenter and faceAreaNormal
    faceAreaNormals_.resize(faces_.size());
    faceCenters_.resize(faces_.size());
    for (size_t fi = 0; fi < faces_.size(); fi++)
    {
        const Face &face = faces_[fi];
        const int nPoints = int(face.size());
        if (nPoints == 3)
        {
            faceCenters_[fi] = triCenter(nodeCoord_[face[0]], nodeCoord_[face[1]], nodeCoord_[face[2]]);
            faceAreaNormals_[fi] = triAreaNormal(nodeCoord_[face[0]], nodeCoord_[face[1]], nodeCoord_[face[2]]);
        }
        else
        {
            Point fCenter = nodeCoord_[face[0]];
            for (int i = 1; i < nPoints; i++)
            {
                fCenter += nodeCoord_[face[i]];
            }
            fCenter /= nPoints;

            Vec sumN(0.0, 0.0, 0.0);
            Vec sumAc(0.0, 0.0, 0.0);
            double sumA = 0.0;

            for (int i = 0; i < nPoints; i++)
            {
                const Point thisPoint = nodeCoord_[face[i]];
                const Point nextPoint = nodeCoord_[i == (nPoints - 1) ? face[0] : face[i + 1]];
                Vec c = thisPoint + nextPoint + fCenter;
                Vec n = (nextPoint - thisPoint) ^ (fCenter - thisPoint);
                double a = mag(n);

                sumN += n;
                sumA += a;
                sumAc += a * c;
            }

            if (sumA < 1e-18)
            {
                faceCenters_[fi] = fCenter;
                faceAreaNormals_[fi] = Vec(0.0, 0.0, 0.0);
            }
            else
            {
                faceCenters_[fi] = (1.0 / 3.0) * sumAc / sumA;
                faceAreaNormals_[fi] = 0.5 * sumN;
            }
        }
    } // end for face data

    // calc cellCenter cellVolume and cellFaces
    PointList geoCellCenter(nCell_, Point(0.0, 0.0, 0.0));
    cellFaces_.resize(nCell_);
    for (size_t i = 0; i < faces_.size(); i++)
    {
        const int nighebor1 = faceNeighbors_[2 * i];
        const int nighebor2 = faceNeighbors_[2 * i + 1];
        if (nighebor1 != -1)
        {
            geoCellCenter[nighebor1] += faceCenters_[i];
            cellFaces_[nighebor1].push_back(int(i));
        }
        if (nighebor2 != -1)
        {
            geoCellCenter[nighebor2] += faceCenters_[i];
            cellFaces_[nighebor2].push_back(int(i));
        }
    }

    for (size_t i = 0; i < geoCellCenter.size(); i++)
    {
        geoCellCenter[i] /= cellFaces_[i].size();
    }

    if (cellCenterType_ == CellCenterType::GeometricCenter)
    {
        cellCenters_ = geoCellCenter;
    }
    else if (cellCenterType_ == CellCenterType::Barycenter)
    {   
        // for baryCellCenter
        PointList baryCellCenter(nCell_, Point(0.0, 0.0, 0.0));
        cellVolumes_.resize(nCell_, 0.0);
        for (size_t i = 0; i < faces_.size(); i++)
        {
            const Point fC(faceCenters_[i]);
            const Vec fA(faceAreaNormals_[i]);
            const int neighbor1 = faceNeighbors_[2 * i];
            const int neighbor2 = faceNeighbors_[2 * i + 1];
            if (neighbor1 != -1)
            {
                const double pyr3Vol = abs(fA & (fC - geoCellCenter[neighbor1]));
                const Vec pc = (3.0 / 4.0) * fC + (1.0 / 4.0) * geoCellCenter[neighbor1];
                baryCellCenter[neighbor1] += pyr3Vol * pc;
                cellVolumes_[neighbor1] += pyr3Vol;
            }
            if (neighbor2 != -1)
            {
                const double pyr3Vol = abs(fA & (fC - geoCellCenter[neighbor2]));
                const Vec pc = (3.0 / 4.0) * fC + (1.0 / 4.0) * geoCellCenter[neighbor2];
                baryCellCenter[neighbor2] += pyr3Vol * pc;
                cellVolumes_[neighbor2] += pyr3Vol;
            }
        }
        cellCenters_.resize(nCell_);
        for (size_t i = 0; i < baryCellCenter.size(); i++)
        {   
            // avoid divide by zero
            if (abs(cellVolumes_[i]) > 1e-20)
            {
                cellCenters_[i] = (baryCellCenter[i] / cellVolumes_[i]);
            }
            else
            {
                cellCenters_[i] = geoCellCenter[i];
            }
            cellVolumes_[i] *= (1.0 / 3.0);
        }
    }

    // calc cellNodes
    cellNodes_.resize(nCell_, IntList());
    for (size_t i = 0; i < cellNodes_.size(); i++)
    {

        IntSet cellNodeSet;
        for (size_t j = 0; j < cellFaces_[i].size(); j++)
        {
            const Face &face = faces_[cellFaces_[i][j]];
            for (size_t k = 0; k < face.size(); k++)
            {
                cellNodeSet.insert(face[k]);
            }
        }
        cellNodes_[i].insert(cellNodes_[i].end(), cellNodeSet.begin(), cellNodeSet.end());
    }
    // calc cellKNodes
    cellKNodes_.resize(nCell_, IntList(8));
    for (size_t i = 0; i < cellKNodes_.size(); i++)
    {
        Face faceK1;
        Face faceK2;
        IntList kFaceIdx;
        // get k1 and k2 face
        for(size_t fi = 0 ; fi<cellFaces_[i].size();fi++)
        {   
            // get k1 and k2 face
            if(faceTags_[cellFaces_[i][fi]] == FaceType::KFace)
            {
                kFaceIdx.push_back(cellFaces_[i][fi]);
            }
            
        }
        if(kFaceIdx.size() == 2)
        {   
            // get k1 and k2 face
            faceK1 = faces_[kFaceIdx[0]];
            faceK2 = faces_[kFaceIdx[1]];
        }
        else
        {
            std::cerr << "Error: Not find two KFaces for a cell!" << std::endl;
            std::abort();
        }

        if (faceK1.size() == 4 && faceK2.size() == 4)
        {
            // cellNodes_ 0-3 is k1 face nodes, 4-7 is k2 face nodes
            cellKNodes_[i][0] = faceK1[0];
            cellKNodes_[i][1] = faceK1[1];
            cellKNodes_[i][2] = faceK1[2];
            cellKNodes_[i][3] = faceK1[3];
            cellKNodes_[i][4] = faceK2[0];
            cellKNodes_[i][5] = faceK2[1];
            cellKNodes_[i][6] = faceK2[2];
            cellKNodes_[i][7] = faceK2[3];
        }
        else
        {
            std::cerr << "Error: cellNodes_ size not match!" << std::endl;
            std::abort();
        }
    }
    // calc boundingbox
    boundingBox_.first = Point(1e20, 1e20, 1e20);
    boundingBox_.second = Point(-1e20, -1e20, -1e20);
    for (size_t i = 0; i < nodeCoord_.size(); i++)
    {
        boundingBox_.first.x = std::min(boundingBox_.first.x, nodeCoord_[i].x);
        boundingBox_.first.y = std::min(boundingBox_.first.y, nodeCoord_[i].y);
        boundingBox_.first.z = std::min(boundingBox_.first.z, nodeCoord_[i].z);
        boundingBox_.second.x = std::max(boundingBox_.second.x, nodeCoord_[i].x);
        boundingBox_.second.y = std::max(boundingBox_.second.y, nodeCoord_[i].y);
        boundingBox_.second.z = std::max(boundingBox_.second.z, nodeCoord_[i].z);
    }
    std::cout<<"X Min: "<<boundingBox_.first.x<< " X Max: "<<boundingBox_.second.x<<std::endl;
    std::cout<<"Y Min: "<<boundingBox_.first.y<< " Y Max: "<<boundingBox_.second.y<<std::endl;
    std::cout<<"Z Min: "<<boundingBox_.first.z<< " Z Max: "<<boundingBox_.second.z<<std::endl;

    // get non-slip fault face
    getCommNormFaultFace();
    
    clock_t endT1 = std::clock();
    if (perfFlag_)
    {
        std::cout << "calcPolyGridData Time: " << double(endT1 - startT) / CLOCKS_PER_SEC << " s" << std::endl;
    }
}
void GRDECLGrid::getNonSlipFaultFace()
{
    // find unregular face
    //Map<int,int> tetFaceAndNodeMap;  //<faceIdx,nodeIdx>
    //Map<int, IntList> tetCellAndFace; //<cellIdx,faceIdxList>
    IntSet addedFace;
    for(size_t i = 0; i < faultFaces_.size(); i++)
    {
        addedFace.insert(faultFaces_[i]);
    }

    int nonQuadFace = 0;

    for (size_t i = 0; i < faces_.size(); i++)
    {
        const Face &face = faces_[i];
        int parallelEdge = 0;
        if (face.size() != 4)
        {
            if (faceTags_[i] == FaceType::KFace)
            {
                std::cerr << "Error grids: k face is not a quad!" << std::endl;
                std::abort();
            }
            if(addedFace.find(int(i)) == addedFace.end())
            {
                addedFace.insert(int(i));
                faultFaces_.push_back(int(i));
            }
            nonQuadFace++;
        }
        for (size_t j = 0; j < face.size(); j++)
        {
            const int nodeIdxPre = face[j == 0 ? face.size() - 1 : j - 1];
            const int nodeIdxCurr = face[j];
            const int nodeIdxNext = face[j == (face.size() - 1) ? 0 : j + 1];
            const Vec v1 = nodeCoord_[nodeIdxCurr] - nodeCoord_[nodeIdxPre];
            const Vec v2 = nodeCoord_[nodeIdxNext] - nodeCoord_[nodeIdxCurr];
            if (coLineVec(v1, v2))
            {
                parallelEdge++;
                // find cell neighbor
                const int cellIdx = faceNeighbors_[2 * i] != -1 ? faceNeighbors_[2 * i] : faceNeighbors_[2 * i + 1];
                if(cellIdx == -1)
                {
                    std::cerr << "Error faces: fault face has only one cell neighbor!" << std::endl;
                    std::abort();
                }
                IntList& cellIFaces = cellFaces_[cellIdx];
                if(cellIFaces.size() != 6)
                {
                    std::cerr << "Error grids: cell do not has 6 faces!" << std::endl;
                    //std::abort();
                }
                // 4 condidtions for non-slip fault face
                if (j == 1)
                {
                    // 0,2 face is non-slip fault face
                    if (addedFace.find(cellIFaces[0]) == addedFace.end())
                    {
                        addedFace.insert(cellIFaces[0]);
                        faultFaces_.push_back(cellIFaces[0]);
                    }
                    if (addedFace.find(cellIFaces[2]) == addedFace.end())
                    {
                        addedFace.insert(cellIFaces[2]);
                        faultFaces_.push_back(cellIFaces[2]);
                    }
                }
                else if (j == 2)
                {
                    // 0,3 face is non-slip fault face
                    if (addedFace.find(cellIFaces[0]) == addedFace.end())
                    {
                        addedFace.insert(cellIFaces[0]);
                        faultFaces_.push_back(cellIFaces[0]);
                    }
                    if (addedFace.find(cellIFaces[3]) == addedFace.end())
                    {
                        addedFace.insert(cellIFaces[3]);
                        faultFaces_.push_back(cellIFaces[3]);
                    }

                }
                else if (j == 3)
                {
                    // 1,3 face is non-slip fault face
                    if (addedFace.find(cellIFaces[1]) == addedFace.end())
                    {
                        addedFace.insert(cellIFaces[1]);
                        faultFaces_.push_back(cellIFaces[1]);
                    }
                    if (addedFace.find(cellIFaces[3]) == addedFace.end())
                    {
                        addedFace.insert(cellIFaces[3]);
                        faultFaces_.push_back(cellIFaces[3]);
                    }
                }
                else if (j == 0)
                {
                    // 1,2 face is non-slip fault face
                    if (addedFace.find(cellIFaces[1]) == addedFace.end())
                    {
                        addedFace.insert(cellIFaces[1]);
                        faultFaces_.push_back(cellIFaces[1]);
                    }
                    if (addedFace.find(cellIFaces[2]) == addedFace.end())
                    {
                        addedFace.insert(cellIFaces[2]);
                        faultFaces_.push_back(cellIFaces[2]);
                    }
                }
            }
        }
        if (parallelEdge > 1)
        {
            std::cerr << "Error grids: face has more than one parallel edge pairs!" << std::endl;
            std::abort();
        }
    }

}

void GRDECLGrid::getCommNormFaultFace()
{
    // find unregular face
    //Map<int,int> tetFaceAndNodeMap;  //<faceIdx,nodeIdx>
    //Map<int, IntList> tetCellAndFace; //<cellIdx,faceIdxList>
    IntSet addedFace;
    // for(size_t i = 0; i < faultFaces_.size(); i++)
    // {
    //     addedFace.insert(faultFaces_[i]);
    // }
    faultFaces_.resize(0);

    int nonQuadFace = 0;

    for (int ci = 0; ci < cellFaces_.size(); ci++)
    {
        for (int f1 = 0; f1 < cellFaces_[ci].size(); f1++)
        {
            IntList &face1 = faces_[cellFaces_[ci][f1]];
            for (int f2 = f1 + 1; f2 < cellFaces_[ci].size(); f2++)
            {
                IntList &face2 = faces_[cellFaces_[ci][f2]];
                if (coEdgeFace(face1, face2) && 
                coPlaneFace(faceAreaNormals_[cellFaces_[ci][f1]], faceAreaNormals_[cellFaces_[ci][f2]]))
                {
                    if (addedFace.find(cellFaces_[ci][f1]) == addedFace.end())
                    {
                        addedFace.insert(cellFaces_[ci][f1]);
                        faultFaces_.push_back(cellFaces_[ci][f1]);
                    }
                    if (addedFace.find(cellFaces_[ci][f2]) == addedFace.end())
                    {
                        addedFace.insert(cellFaces_[ci][f2]);
                        faultFaces_.push_back(cellFaces_[ci][f2]);
                    }
                }
            }
        }
    }
}

void GRDECLGrid::decomposeAllTetCell()
{   
    IMinFaceNodes_.clear();
    IMaxFaceNodes_.clear();
    JMinFaceNodes_.clear();
    JMaxFaceNodes_.clear();
    KMinFaceNodes_.clear();
    KMaxFaceNodes_.clear();
    clock_t startT = std::clock();
    allTetCoord_ = nodeCoord_;
    // find unregular K face
    Map<int, int> tetFace;//<faceIdx,nodeIdx>
    IntSet tetFaceCell; //<cellIdx,faceIdxList>
    for (size_t i = 0; i < faces_.size(); i++)
    {   
        const Face &face = faces_[i];
        // if (face.size() != 4)
        // {
        //     std::cerr << "Error grids: face is not quad!" << std::endl;
        // }
        //         // skip non-K face
        if(faceTags_[i] != FaceType::KFace)
        {
            continue;
        }
        int parallelEdge = 0;
        
        for (size_t j = 0; j < face.size(); j++)
        {
            // find pinched face
            const int nodeIdxPre = face[j == 0 ? face.size() - 1 : j - 1];
            const int nodeIdxCurr = face[j];
            const int nodeIdxNext = face[j == (face.size() - 1) ? 0 : j + 1];
            const Vec v1 = allTetCoord_[nodeIdxCurr] - allTetCoord_[nodeIdxPre];
            const Vec v2 = allTetCoord_[nodeIdxNext] - allTetCoord_[nodeIdxCurr];
            if (coLineVec(v1, v2))
            {
                parallelEdge++;
                // pinched face
                if (j == 0 || j == 2)
                {
                    tetFace.emplace(int(i), int(j));
                    if (faceNeighbors_[2 * i] != -1)
                    {
                        tetFaceCell.insert(faceNeighbors_[2 * i]);
                        IntList &cellFace = cellFaces_[faceNeighbors_[2 * i]];
                        int anotherKFace = -1;
                        for (size_t k = 0; k < cellFace.size(); k++)
                        {
                            if (faceTags_[cellFace[k]] == FaceType::KFace && cellFace[k] != i)
                            {
                                anotherKFace = cellFace[k];
                                break;
                            }
                        }
                        if (anotherKFace == -1)
                        {
                            std::cerr << "Error grids: cell do not has another K face!" << std::endl;
                            std::abort();
                        }
                        tetFace.emplace(anotherKFace, int(j));
                    }
                    if (faceNeighbors_[2 * i + 1] != -1)
                    {
                        tetFaceCell.insert(faceNeighbors_[2 * i + 1]);
                        IntList &cellFace = cellFaces_[faceNeighbors_[2 * i + 1]];
                        int anotherKFace = -1;
                        for (size_t k = 0; k < cellFace.size(); k++)
                        {
                            if (faceTags_[cellFace[k]] == FaceType::KFace && cellFace[k] != i)
                            {
                                anotherKFace = cellFace[k];
                                break;
                            }
                        }
                        if (anotherKFace == -1)
                        {
                            std::cerr << "Error grids: cell do not has another K face!" << std::endl;
                            std::abort();
                        }
                        tetFace.emplace(anotherKFace, int(j));
                    }
                }
            }
        }
        if (parallelEdge > 1)
        {
            std::cerr << "Error grids: face has more than one parallel edge pairs!" << std::endl;
            std::abort();
        }
    }

    // select boundary face
    for (int i = 0; i < int(faces_.size()); i++)
    {   
        
        if (IMinFace_.find(i) != IMinFace_.end())
        {
            IMinFaceNodes_.push_back(faces_[i]);
            // IMinFaceNodes_.push_back(IntList{faces_[i][1],faces_[i][3],faces_[i][2]});
            // IMinFaceNodes_.push_back(IntList{faces_[i][0],faces_[i][1],faces_[i][3]});
        }
        else if (IMaxFace_.find(i) != IMaxFace_.end())
        {
            IMaxFaceNodes_.push_back(faces_[i]);
            // IMaxFaceNodes_.push_back(IntList{faces_[i][1],faces_[i][3],faces_[i][2]});
            // IMaxFaceNodes_.push_back(IntList{faces_[i][0],faces_[i][1],faces_[i][3]});
        }
        else if (JMinFace_.find(i) != JMinFace_.end())
        {
            JMinFaceNodes_.push_back(faces_[i]);
            // JMinFaceNodes_.push_back(IntList{faces_[i][1],faces_[i][3],faces_[i][2]});
            // JMinFaceNodes_.push_back(IntList{faces_[i][0],faces_[i][1],faces_[i][3]});
        }
        else if (JMaxFace_.find(i) != JMaxFace_.end())
        {
            JMaxFaceNodes_.push_back(faces_[i]);
            // JMaxFaceNodes_.push_back(IntList{faces_[i][1],faces_[i][3],faces_[i][2]});
            // JMaxFaceNodes_.push_back(IntList{faces_[i][0],faces_[i][1],faces_[i][3]});
        }
        else if (KMinFace_.find(i) != KMinFace_.end())
        {   
            KMinFaceNodes_.push_back(faces_[i]);
            // if(tetFace.find(i) == tetFace.end())
            // {
            //     // KMinFaceNodes_.push_back(IntList{faces_[i][1],faces_[i][3],faces_[i][2]});
            //     // KMinFaceNodes_.push_back(IntList{faces_[i][0],faces_[i][1],faces_[i][3]});
            // }
            // else
            // {
            //     // KMinFaceNodes_.push_back(IntList{faces_[i][0],faces_[i][1],faces_[i][2]});
            //     // KMinFaceNodes_.push_back(IntList{faces_[i][0],faces_[i][2],faces_[i][3]});  
            // }
        }
        else if (KMaxFace_.find(i) != KMaxFace_.end())
        {
            KMaxFaceNodes_.push_back(faces_[i]);
            // if(tetFace.find(i) == tetFace.end())
            // {
            //     // KMaxFaceNodes_.push_back(IntList{faces_[i][1],faces_[i][3],faces_[i][2]});
            //     // KMaxFaceNodes_.push_back(IntList{faces_[i][0],faces_[i][1],faces_[i][3]});
            // }
            // else
            // {
            //     // KMaxFaceNodes_.push_back(IntList{faces_[i][0],faces_[i][1],faces_[i][2]});
            //     // KMaxFaceNodes_.push_back(IntList{faces_[i][0],faces_[i][2],faces_[i][3]});                  
            // }
        }
    }

    for (int hexI = 0; hexI < nCell_; hexI++)
    {
        const IntList &hexNodes = cellKNodes_[hexI];
        // standard decompose
        if (tetFaceCell.find(hexI) == tetFaceCell.end())
        {
            // half based on 0 node
            allTetCellNodes_.emplace_back(IntList{hexNodes[0], hexNodes[1], hexNodes[3], hexNodes[5]});
            allTetCellNodes_.emplace_back(IntList{hexNodes[0], hexNodes[3], hexNodes[5], hexNodes[7]});
            allTetCellNodes_.emplace_back(IntList{hexNodes[0], hexNodes[4], hexNodes[7], hexNodes[5]});
            // half based on 6 node
            allTetCellNodes_.emplace_back(IntList{hexNodes[6], hexNodes[2], hexNodes[1], hexNodes[3]});
            allTetCellNodes_.emplace_back(IntList{hexNodes[6], hexNodes[5], hexNodes[3], hexNodes[1]});
            allTetCellNodes_.emplace_back(IntList{hexNodes[6], hexNodes[7], hexNodes[3], hexNodes[5]});
        }
        else // alternative decompose
        {
            //  decompose prism cell to 6 tet cells
            allTetCellNodes_.emplace_back(IntList{hexNodes[0], hexNodes[1], hexNodes[2], hexNodes[6]});
            allTetCellNodes_.emplace_back(IntList{hexNodes[0], hexNodes[1], hexNodes[5], hexNodes[6]});
            allTetCellNodes_.emplace_back(IntList{hexNodes[0], hexNodes[4], hexNodes[5], hexNodes[6]});
            //
            allTetCellNodes_.emplace_back(IntList{hexNodes[0], hexNodes[2], hexNodes[3], hexNodes[6]});
            allTetCellNodes_.emplace_back(IntList{hexNodes[0], hexNodes[3], hexNodes[7], hexNodes[6]});
            allTetCellNodes_.emplace_back(IntList{hexNodes[0], hexNodes[4], hexNodes[7], hexNodes[6]});
        }
    }

    tetReserviorCellNum_ = int(allTetCellNodes_.size());
    clock_t endT = std::clock();
    if (perfFlag_)
    {
        std::cout << "Cell TetDecompose Time: " << double(endT - startT) / CLOCKS_PER_SEC << " s" << std::endl;
    }
}

void GRDECLGrid::decomposeAllTetExtendCell(){
    // data :: allTetCoord_,allTetCellNodes_,tetPointData_,tetCellData_,tetReserviorCellNum_
    IntListList allTetExtendCellNodes;
    IntList tetExtendCellData;
    for (int i = tetReserviorCellNum_; i < int(allTetCellNodes_.size()); i++)
    {
        allTetExtendCellNodes.push_back(allTetCellNodes_[i]);
        tetExtendCellData.push_back(tetCellData_[i]);
    }
    allTetCellNodes_.resize(tetReserviorCellNum_);
    tetCellData_.resize(tetReserviorCellNum_);
    const array<int, 8> rIdx = {4,5,6,7,0,1,2,3};
    for(int i =0 ;i<allTetExtendCellNodes.size();i++)
    {
        IntList &hexNodes = allTetExtendCellNodes[i];
        if (tetExtendCellData[i] == 3 || tetExtendCellData[i] == 5)
        {
            Face kFace{hexNodes[0], hexNodes[1], hexNodes[2], hexNodes[3]};
            bool altDecomp = false;
            for (size_t j = 0; j < kFace.size(); j++)
            {
                // find pinched point
                const int nodeIdxPre = kFace[j == 0 ? kFace.size() - 1 : j - 1];
                const int nodeIdxCurr = kFace[j];
                const int nodeIdxNext = kFace[j == (kFace.size() - 1) ? 0 : j + 1];
                const Vec v1 = allTetCoord_[nodeIdxCurr] - allTetCoord_[nodeIdxPre];
                const Vec v2 = allTetCoord_[nodeIdxNext] - allTetCoord_[nodeIdxCurr];
                const Vec v1Norm = v1 / mag(v1);
                const Vec v2Norm = v2 / mag(v2);
                if (coLineVec(v1, v2))
                {
                    // pinched face
                    if (j == 0 || j == 2)
                    {
                        altDecomp = true;
                    }
                }
            }
            if (altDecomp)
            {
                if (tetExtendCellData[i] == 3)
                {
                    //  decompose prism cell to 6 tet cells
                    allTetCellNodes_.emplace_back(IntList{hexNodes[0], hexNodes[1], hexNodes[2], hexNodes[6]});
                    allTetCellNodes_.emplace_back(IntList{hexNodes[0], hexNodes[1], hexNodes[5], hexNodes[6]});
                    allTetCellNodes_.emplace_back(IntList{hexNodes[0], hexNodes[4], hexNodes[5], hexNodes[6]});
                    //
                    allTetCellNodes_.emplace_back(IntList{hexNodes[0], hexNodes[2], hexNodes[3], hexNodes[6]});
                    allTetCellNodes_.emplace_back(IntList{hexNodes[0], hexNodes[3], hexNodes[7], hexNodes[6]});
                    allTetCellNodes_.emplace_back(IntList{hexNodes[0], hexNodes[4], hexNodes[7], hexNodes[6]});
                }
                else if (tetExtendCellData[i] == 5)
                {
                    allTetCellNodes_.emplace_back(IntList{hexNodes[rIdx[0]], hexNodes[rIdx[1]], hexNodes[rIdx[2]], hexNodes[rIdx[6]]});
                    allTetCellNodes_.emplace_back(IntList{hexNodes[rIdx[0]], hexNodes[rIdx[1]], hexNodes[rIdx[5]], hexNodes[rIdx[6]]});
                    allTetCellNodes_.emplace_back(IntList{hexNodes[rIdx[0]], hexNodes[rIdx[4]], hexNodes[rIdx[5]], hexNodes[rIdx[6]]});
                    //
                    allTetCellNodes_.emplace_back(IntList{hexNodes[rIdx[0]], hexNodes[rIdx[2]], hexNodes[rIdx[3]], hexNodes[rIdx[6]]});
                    allTetCellNodes_.emplace_back(IntList{hexNodes[rIdx[0]], hexNodes[rIdx[3]], hexNodes[rIdx[7]], hexNodes[rIdx[6]]});
                    allTetCellNodes_.emplace_back(IntList{hexNodes[rIdx[0]], hexNodes[rIdx[4]], hexNodes[rIdx[7]], hexNodes[rIdx[6]]});
                }
            tetCellData_.insert(tetCellData_.end(), 6, tetExtendCellData[i]);
            }
            else
            {
                if (tetExtendCellData[i] == 3)
                {
                    // half based on 0 node
                    allTetCellNodes_.emplace_back(IntList{hexNodes[0], hexNodes[1], hexNodes[3], hexNodes[5]});
                    allTetCellNodes_.emplace_back(IntList{hexNodes[0], hexNodes[3], hexNodes[5], hexNodes[7]});
                    allTetCellNodes_.emplace_back(IntList{hexNodes[0], hexNodes[4], hexNodes[7], hexNodes[5]});
                    // half based on 6 node
                    allTetCellNodes_.emplace_back(IntList{hexNodes[6], hexNodes[2], hexNodes[1], hexNodes[3]});
                    allTetCellNodes_.emplace_back(IntList{hexNodes[6], hexNodes[5], hexNodes[3], hexNodes[1]});
                    allTetCellNodes_.emplace_back(IntList{hexNodes[6], hexNodes[7], hexNodes[3], hexNodes[5]});
                }
                else if (tetExtendCellData[i] == 5)
                {
                    // half based on 0 node
                    allTetCellNodes_.emplace_back(IntList{hexNodes[rIdx[0]], hexNodes[rIdx[1]], hexNodes[rIdx[3]], hexNodes[rIdx[5]]});
                    allTetCellNodes_.emplace_back(IntList{hexNodes[rIdx[0]], hexNodes[rIdx[3]], hexNodes[rIdx[5]], hexNodes[rIdx[7]]});
                    allTetCellNodes_.emplace_back(IntList{hexNodes[rIdx[0]], hexNodes[rIdx[4]], hexNodes[rIdx[7]], hexNodes[rIdx[5]]});
                    // half based on 6 node
                    allTetCellNodes_.emplace_back(IntList{hexNodes[rIdx[6]], hexNodes[rIdx[2]], hexNodes[rIdx[1]], hexNodes[rIdx[3]]});
                    allTetCellNodes_.emplace_back(IntList{hexNodes[rIdx[6]], hexNodes[rIdx[5]], hexNodes[rIdx[3]], hexNodes[rIdx[1]]});
                    allTetCellNodes_.emplace_back(IntList{hexNodes[rIdx[6]], hexNodes[rIdx[7]], hexNodes[rIdx[3]], hexNodes[rIdx[5]]});
                }
                tetCellData_.insert(tetCellData_.end(), 6, tetExtendCellData[i]);
            }
        }
        else
        {
            //reorderCellNodes(hexNodes, cIdx);
            // half based on 0 node
            // allTetCellNodes_.emplace_back(IntList{hexNodes[cIdx[0]], hexNodes[cIdx[1]], hexNodes[cIdx[3]], hexNodes[cIdx[5]]});
            // allTetCellNodes_.emplace_back(IntList{hexNodes[cIdx[0]], hexNodes[cIdx[3]], hexNodes[cIdx[5]], hexNodes[cIdx[7]]});
            // allTetCellNodes_.emplace_back(IntList{hexNodes[cIdx[0]], hexNodes[cIdx[4]], hexNodes[cIdx[7]], hexNodes[cIdx[5]]});
            // // half based on 6 node
            // allTetCellNodes_.emplace_back(IntList{hexNodes[cIdx[6]], hexNodes[cIdx[2]], hexNodes[cIdx[1]], hexNodes[cIdx[3]]});
            // allTetCellNodes_.emplace_back(IntList{hexNodes[cIdx[6]], hexNodes[cIdx[5]], hexNodes[cIdx[3]], hexNodes[cIdx[1]]});
            // allTetCellNodes_.emplace_back(IntList{hexNodes[cIdx[6]], hexNodes[cIdx[7]], hexNodes[cIdx[3]], hexNodes[cIdx[5]]});
            reorderCellNodes(allTetCoord_, hexNodes);
            allTetCellNodes_.emplace_back(IntList{hexNodes[0], hexNodes[1], hexNodes[3], hexNodes[5]});
            allTetCellNodes_.emplace_back(IntList{hexNodes[0], hexNodes[3], hexNodes[5], hexNodes[7]});
            allTetCellNodes_.emplace_back(IntList{hexNodes[0], hexNodes[4], hexNodes[7], hexNodes[5]});
            // half based on 6 node
            allTetCellNodes_.emplace_back(IntList{hexNodes[6], hexNodes[2], hexNodes[1], hexNodes[3]});
            allTetCellNodes_.emplace_back(IntList{hexNodes[6], hexNodes[5], hexNodes[3], hexNodes[1]});
            allTetCellNodes_.emplace_back(IntList{hexNodes[6], hexNodes[7], hexNodes[3], hexNodes[5]});
            tetCellData_.insert(tetCellData_.end(), 6, tetExtendCellData[i]);
        }
    }
}

void GRDECLGrid::decomposeHybridExtendCell(){
    for (int i = decompReservoirCellNum_; i < decomposedCellNodes_.size(); i++)
    {
        // decompose prism cell to 3 pyramid cells and 2 tet cells
        if (decomposedCellNodes_[i].size() == 6)
        {   
            
            const int centerIdx = int(decomposedCoord_.size());
            IntList prismNodes = decomposedCellNodes_[i];
            decomposedCoord_.push_back(PointsCenter(decomposedCoord_,prismNodes));
            decomposedPointData_.push_back(0);
            const int prismCellData = decomposedCellData_[i];
            // pyramid 1
            // for two tet
            decomposedCellNodes_[i] = IntList{prismNodes[0], prismNodes[1], prismNodes[2], centerIdx};
            decomposedCellNodes_.emplace_back(IntList{prismNodes[3], prismNodes[4], prismNodes[5], centerIdx});
            // for three pyramid
            decomposedCellNodes_.emplace_back(IntList{prismNodes[0], prismNodes[1], prismNodes[4], prismNodes[3], centerIdx});
            decomposedCellNodes_.emplace_back(IntList{prismNodes[1], prismNodes[2], prismNodes[5], prismNodes[4], centerIdx});
            decomposedCellNodes_.emplace_back(IntList{prismNodes[0], prismNodes[2], prismNodes[5], prismNodes[3], centerIdx});
            decomposedCellData_.insert(decomposedCellData_.end(), 4, prismCellData);
        }
    }
}

void GRDECLGrid::decomposeAllHexCell()
{   
    IMinFaceNodes_.clear();
    IMaxFaceNodes_.clear();
    JMinFaceNodes_.clear();
    JMaxFaceNodes_.clear();
    KMinFaceNodes_.clear();
    KMaxFaceNodes_.clear();
    allHexCoord_ = nodeCoord_;
    clock_t startT = std::clock();
    // find unregular face
    Map<int,int> tetFaceAndNodeMap;  //<faceIdx,nodeIdx>
    vector<IntPair> tetFaceAndNode;   //<faceIdx,nodeIdx>
    Map<int, IntList> tetCellAndFace; //<cellIdx,faceIdxList>
    for (size_t i = 0; i < faces_.size(); i++)
    {
        const Face &face = faces_[i];
        int parallelEdge = 0;
        if(tetFaceAndNodeMap.find(int(i)) != tetFaceAndNodeMap.end())
        {
            continue;
        }
        if(faceTags_[i] != FaceType::KFace)
        {
            continue;
        }
        for (size_t j = 0; j < face.size(); j++)
        {   
            // find pinched face
            const int nodeIdxPre = face[j == 0 ? face.size() - 1 : j - 1];
            const int nodeIdxCurr = face[j];
            const int nodeIdxNext = face[j == (face.size() - 1) ? 0 : j + 1];
            const Vec v1 = allHexCoord_[nodeIdxCurr] - allHexCoord_[nodeIdxPre];
            const Vec v2 = allHexCoord_[nodeIdxNext] - allHexCoord_[nodeIdxCurr];
            if (coLineVec(v1, v2))
            {
                parallelEdge++;
                tetFaceAndNode.push_back({int(i), int(j)});
                tetFaceAndNodeMap.emplace(int(i), int(j));
                if (faceNeighbors_[2 * i] != -1)
                {
                    IntList& cellFace = cellFaces_[faceNeighbors_[2 * i]];
                    int anotherKFace = -1;
                    for(size_t k = 0; k < cellFace.size(); k++)
                    {
                        if(faceTags_[cellFace[k]] == FaceType::KFace && cellFace[k] != i)
                        {
                            anotherKFace = cellFace[k];
                            break;
                        }
                    }
                    if(anotherKFace == -1)
                    {
                        std::cerr << "Error grids: cell do not has another K face!" << std::endl;
                        std::abort();
                    }
                    if (tetCellAndFace.find(faceNeighbors_[2 * i]) == tetCellAndFace.end())
                    {
                        tetCellAndFace.emplace(faceNeighbors_[2 * i], IntList());
                    }
                    tetCellAndFace[faceNeighbors_[2 * i]].push_back(int(i));
                    tetCellAndFace[faceNeighbors_[2 * i]].push_back(anotherKFace);
                    if(tetCellAndFace[faceNeighbors_[2 * i]][0]>tetCellAndFace[faceNeighbors_[2 * i]][1])
                    {
                        std::swap(tetCellAndFace[faceNeighbors_[2 * i]][0],tetCellAndFace[faceNeighbors_[2 * i]][1]);
                    }
                    tetFaceAndNode.push_back({anotherKFace, int(j)});
                    tetFaceAndNodeMap.emplace(anotherKFace, int(j));
                }
                if (faceNeighbors_[2 * i + 1] != -1)
                {   
                    IntList& cellFace = cellFaces_[faceNeighbors_[2 * i + 1]];
                    int anotherKFace = -1;
                    for(size_t k = 0; k < cellFace.size(); k++)
                    {
                        if(faceTags_[cellFace[k]] == FaceType::KFace && cellFace[k] != i)
                        {
                            anotherKFace = cellFace[k];
                            break;
                        }
                    }
                    if(anotherKFace == -1)
                    {
                        std::cerr << "Error grids: cell do not has another K face!" << std::endl;
                        std::abort();
                    }
                    if (tetCellAndFace.find(faceNeighbors_[2 * i + 1]) == tetCellAndFace.end())
                    {
                        tetCellAndFace.emplace(faceNeighbors_[2 * i + 1], IntList());
                    }
                    tetCellAndFace[faceNeighbors_[2 * i + 1]].push_back(int(i));
                    tetCellAndFace[faceNeighbors_[2 * i + 1]].push_back(anotherKFace);
                    if (tetCellAndFace[faceNeighbors_[2 * i + 1]][0] > tetCellAndFace[faceNeighbors_[2 * i + 1]][1])
                    {
                        std::swap(tetCellAndFace[faceNeighbors_[2 * i + 1]][0], tetCellAndFace[faceNeighbors_[2 * i + 1]][1]);
                    }
                    tetFaceAndNode.push_back({anotherKFace, int(j)});
                    tetFaceAndNodeMap.emplace(anotherKFace, int(j));
                }
            }
        }
        if (parallelEdge > 1)
        {
            std::cerr << "Error grids: face has more than one parallel edge pairs!" << std::endl;
            std::abort();
        }
    }
    // select boundary face
    for (int i = 0; i < int(faces_.size()); i++)
    {

        if (IMinFace_.find(i) != IMinFace_.end())
        {
            IMinFaceNodes_.push_back(faces_[i]);
        }
        else if (IMaxFace_.find(i) != IMaxFace_.end())
        {
            IMaxFaceNodes_.push_back(faces_[i]);
        }
        else if (JMinFace_.find(i) != JMinFace_.end())
        {
            JMinFaceNodes_.push_back(faces_[i]);
        }
        else if (JMaxFace_.find(i) != JMaxFace_.end())
        {
            JMaxFaceNodes_.push_back(faces_[i]);
        }
        else if (KMinFace_.find(i) != KMinFace_.end() &&
                 tetFaceAndNodeMap.find(i) == tetFaceAndNodeMap.end())
        {
            KMinFaceNodes_.push_back(faces_[i]);
        }
        else if (KMaxFace_.find(i) != KMaxFace_.end() &&
                 tetFaceAndNodeMap.find(i) == tetFaceAndNodeMap.end())
        {
            KMaxFaceNodes_.push_back(faces_[i]);
        }
    }
    Map<int, array<int, 8>> decomposedFace;
    for (size_t i = 0; i < tetFaceAndNode.size(); i++)
    {   
        // get face and node
        const int faceIdx = tetFaceAndNode[i].first;
        const Face &faceI = faces_[faceIdx];
        const Point faceCenter = faceCenters_[faceIdx];
        const int tetNodeLocalIdx = tetFaceAndNode[i].second;
        const int tetNodeIdx = faceI[tetNodeLocalIdx];
        const Point tetNode = allHexCoord_[tetNodeIdx];
        const int tetNodeNextLocalIdx = nextIdx(tetNodeLocalIdx, int(faceI.size()));
        const int tetNodeNextIdx = faceI[tetNodeNextLocalIdx];
        const Point tetNodeNext = allHexCoord_[tetNodeNextIdx];
        const int tetNodeOppoLocalIdx = nextIdx(tetNodeNextLocalIdx, int(faceI.size()));
        const int tetNodeOppoIdx = faceI[tetNodeOppoLocalIdx];
        const Point tetNodeOppo = allHexCoord_[tetNodeOppoIdx];
        const int tetNodePreLocalIdx = prevIdx(tetNodeLocalIdx, int(faceI.size()));
        const int tetNodePreIdx = faceI[tetNodePreLocalIdx];
        const Point tetNodePre = allHexCoord_[tetNodePreIdx];
        // calculate new node on edge
        const Point addNodeTet = faceCenter + 0.5 * (tetNode - faceCenter);
        const Point addNodeTetNext = faceCenter + 0.2 * (tetNodeNext - faceCenter);
        const Point addNodeTetOppo = faceCenter + 0.2 * (tetNodeOppo - faceCenter);
        const Point addNodeTetPre = faceCenter + 0.2 * (tetNodePre - faceCenter);

        const int tetNodeAddIdx = int(allHexCoord_.size());
        const int tetNodeAddNextIdx = int(allHexCoord_.size()) + 1;
        const int tetNodeAddOppoIdx = int(allHexCoord_.size()) + 2;
        const int tetNodeAddPreIdx = int(allHexCoord_.size()) + 3;

        decomposedFace.emplace(
            faceIdx,
            array<int, 8>{
                tetNodeAddIdx,
                tetNodeIdx,
                tetNodeAddNextIdx,
                tetNodeNextIdx,
                tetNodeAddOppoIdx,
                tetNodeOppoIdx,
                tetNodeAddPreIdx,
                tetNodePreIdx});

        if(KMinFace_.find(faceIdx) != KMinFace_.end())
        {
            KMinFaceNodes_.push_back(IntList{tetNodeAddIdx, tetNodeAddNextIdx, tetNodeAddOppoIdx, tetNodeAddPreIdx});
            KMinFaceNodes_.push_back(IntList{tetNodeAddIdx, tetNodeIdx, tetNodeNextIdx, tetNodeAddNextIdx});
            KMinFaceNodes_.push_back(IntList{tetNodeAddNextIdx, tetNodeNextIdx, tetNodeOppoIdx, tetNodeAddOppoIdx});
            KMinFaceNodes_.push_back(IntList{tetNodeAddOppoIdx, tetNodeOppoIdx, tetNodePreIdx, tetNodeAddPreIdx});
            KMinFaceNodes_.push_back(IntList{tetNodeAddPreIdx, tetNodePreIdx, tetNodeIdx, tetNodeAddIdx});
        }
        else if(KMaxFace_.find(faceIdx) != KMaxFace_.end())
        {
            KMaxFaceNodes_.push_back(IntList{tetNodeAddIdx, tetNodeAddNextIdx, tetNodeAddOppoIdx, tetNodeAddPreIdx});
            KMaxFaceNodes_.push_back(IntList{tetNodeAddIdx, tetNodeIdx, tetNodeNextIdx, tetNodeAddNextIdx});
            KMaxFaceNodes_.push_back(IntList{tetNodeAddNextIdx, tetNodeNextIdx, tetNodeOppoIdx, tetNodeAddOppoIdx});
            KMaxFaceNodes_.push_back(IntList{tetNodeAddOppoIdx, tetNodeOppoIdx, tetNodePreIdx, tetNodeAddPreIdx});
            KMaxFaceNodes_.push_back(IntList{tetNodeAddPreIdx, tetNodePreIdx, tetNodeIdx, tetNodeAddIdx});
        }
        allHexCoord_.emplace_back(addNodeTet);
        allHexCoord_.emplace_back(addNodeTetNext);
        allHexCoord_.emplace_back(addNodeTetOppo);
        allHexCoord_.emplace_back(addNodeTetPre);

    }
    for (int hexI = 0; hexI < nCell_; hexI++)
    {
        
        // not decompose
        if (tetCellAndFace.find(hexI) == tetCellAndFace.end())
        {
            const IntList &hexNodes = cellKNodes_[hexI];
            allHexCellNodes_.emplace_back(hexNodes);
        }
        else // alternative decompose
        {
            //  decompose prism cell to 5 hex cells
            const IntList &tetFacesList = tetCellAndFace[hexI];
            if(tetFacesList.size() != 2 )
            {
                std::cerr << "Error grids: number of tet face is not match!" << std::endl;
                std::abort();
            }
            const int tetFace0Idx = tetFacesList[0];
            const int tetFace1Idx = tetFacesList[1];
            if(tetFaceAndNodeMap[tetFace0Idx] != tetFaceAndNodeMap[tetFace1Idx])
            {
                std::cerr << "Error grids: two tet face nodeIdx is not match!" << std::endl;
                std::abort();
            }
            const array<int, 8> &face0DecomposedNodes = decomposedFace[tetFacesList[0]];
            const array<int, 8> &face1DecomposedNodes = decomposedFace[tetFacesList[1]];
            // add center hex node
            allHexCellNodes_.emplace_back(
                IntList{face0DecomposedNodes[0], face0DecomposedNodes[2], face0DecomposedNodes[4], face0DecomposedNodes[6],
                        face1DecomposedNodes[0], face1DecomposedNodes[2], face1DecomposedNodes[4], face1DecomposedNodes[6]});
            // add other 4 hexes nodes
            allHexCellNodes_.emplace_back(
                IntList{face0DecomposedNodes[0], face0DecomposedNodes[1], face0DecomposedNodes[3], face0DecomposedNodes[2],
                        face1DecomposedNodes[0], face1DecomposedNodes[1], face1DecomposedNodes[3], face1DecomposedNodes[2]});

            allHexCellNodes_.emplace_back(
                IntList{face0DecomposedNodes[0], face0DecomposedNodes[6], face0DecomposedNodes[7], face0DecomposedNodes[1],
                        face1DecomposedNodes[0], face1DecomposedNodes[6], face1DecomposedNodes[7], face1DecomposedNodes[1]});

            allHexCellNodes_.emplace_back(
                IntList{face0DecomposedNodes[4], face0DecomposedNodes[2], face0DecomposedNodes[3], face0DecomposedNodes[5],
                        face1DecomposedNodes[4], face1DecomposedNodes[2], face1DecomposedNodes[3], face1DecomposedNodes[5]});

            allHexCellNodes_.emplace_back(
                IntList{face0DecomposedNodes[4], face0DecomposedNodes[5], face0DecomposedNodes[7], face0DecomposedNodes[6],
                        face1DecomposedNodes[4], face1DecomposedNodes[5], face1DecomposedNodes[7], face1DecomposedNodes[6]});
        }
    }
    hexReserviorCellNum_ = int(allHexCellNodes_.size());
    clock_t endT = std::clock();
    if (perfFlag_)
    {
        std::cout << "Cell HexDecompose Time: " << double(endT - startT) / CLOCKS_PER_SEC << " s" << std::endl;
    }
}

void GRDECLGrid::decomposeCell()
{   
    IMinFaceNodes_.clear();
    IMaxFaceNodes_.clear();
    JMinFaceNodes_.clear();
    JMaxFaceNodes_.clear();
    KMinFaceNodes_.clear();
    KMaxFaceNodes_.clear();

    clock_t startT = std::clock();
    decomposedCoord_ = nodeCoord_;
    // find unregular face
    IntSet markedCell;
    IntSet markedFace;
    for (size_t i = 0; i < faces_.size(); i++)
    {
        const Face &face = faces_[i];
        // if (face.size() > 4)
        // {
        //     markedFace.insert(int(i));
        //     if (faceNeighbors_[2 * i] != -1)
        //     {
        //         markedCell.insert(faceNeighbors_[2 * i]);
        //     }
        //     if (faceNeighbors_[2 * i + 1] != -1)
        //     {
        //         markedCell.insert(faceNeighbors_[2 * i + 1]);
        //     }
        //     continue;
        // }
        // skip non-K face
        if(faceTags_[i] != FaceType::KFace)
        {
            continue;
        }
        int parallelEdge = 0;
        for (size_t j = 0; j < face.size(); j++)
        {
            // find pinched face
            const int nodeIdxPre = face[j == 0 ? face.size() - 1 : j - 1];
            const int nodeIdxCurr = face[j];
            const int nodeIdxNext = face[j == face.size() - 1 ? 0 : j + 1];
            const Vec v1 = decomposedCoord_[nodeIdxCurr] - decomposedCoord_[nodeIdxPre];
            const Vec v2 = decomposedCoord_[nodeIdxNext] - decomposedCoord_[nodeIdxCurr];
            if (coLineVec(v1, v2))
            {
                parallelEdge++;
                markedFace.insert(int(i));
                if (faceNeighbors_[2 * i] != -1)
                {
                    markedCell.insert(faceNeighbors_[2 * i]);
                }
                if (faceNeighbors_[2 * i + 1] != -1)
                {
                    markedCell.insert(faceNeighbors_[2 * i + 1]);
                }
            }
        }
        if (parallelEdge > int(face.size() - 3))
        {
            std::cerr << "Error grids: face has more than:" << face.size() - 3 << " parallel edge pairs!" << std::endl;
            std::abort();
        }
    }
        // select boundary face
    for (int i = 0; i < int(faces_.size()); i++)
    {

        if (IMinFace_.find(i) != IMinFace_.end())
        {
            IMinFaceNodes_.push_back(faces_[i]);
        }
        else if (IMaxFace_.find(i) != IMaxFace_.end())
        {
            IMaxFaceNodes_.push_back(faces_[i]);
        }
        else if (JMinFace_.find(i) != JMinFace_.end())
        {
            JMinFaceNodes_.push_back(faces_[i]);
        }
        else if (JMaxFace_.find(i) != JMaxFace_.end())
        {
            JMaxFaceNodes_.push_back(faces_[i]);
        }
        else if (KMinFace_.find(i) != KMinFace_.end() &&
                 markedFace.find(i) == markedFace.end())
        {
            KMinFaceNodes_.push_back(faces_[i]);
        }
        else if (KMaxFace_.find(i) != KMaxFace_.end() &&
                 markedFace.find(i) == markedFace.end())
        {
            KMaxFaceNodes_.push_back(faces_[i]);
        }
    }
    Map<int, int> addedFaceAndCenter; //<faceIdx,addedNodeIdx>
    for (size_t ci = 0; ci < cellFaces_.size(); ci++)
    {
        const IntList &facesIdx = cellFaces_[ci];
        // if (facesIdx.size() > 6)
        // {
        //     markedCell.insert(int(ci));
        // }
        if (markedCell.find(int(ci)) != markedCell.end())
        {
            decomposedCoord_.push_back(cellCenters_[ci]);
            const int addedCellCenterIdx = int(decomposedCoord_.size()) - 1;
            for (size_t fi = 0; fi < facesIdx.size(); fi++)
            {   
                
                const int faceID = facesIdx[fi];
                if (faceTags_[faceID] != FaceType::KFace)
                {
                    continue;
                }
                const Face &faceI = faces_[faceID];
                if (markedFace.find(faceID) != markedFace.end())
                { // decompose face
                    IntList newCellNodes(faceI.size() + 1);
                    int addedFaceCenterIdx = -1;
                    if (addedFaceAndCenter.find(faceID) != addedFaceAndCenter.end())
                    {
                        addedFaceCenterIdx = addedFaceAndCenter[faceID];
                    }
                    else
                    {
                        decomposedCoord_.push_back(faceCenters_[faceID]);
                        addedFaceCenterIdx = int(decomposedCoord_.size()) - 1;
                        addedFaceAndCenter.emplace(faceID, addedFaceCenterIdx);
                        for (size_t j = 0; j < faceI.size(); j++)
                        {
                            const int nodeIdxCurr = faceI[j];
                            const int nodeIdxNext = faceI[j == faceI.size() - 1 ? 0 : j + 1];
                            if (KMinFace_.find(faceID) != KMinFace_.end())
                            {
                                KMinFaceNodes_.push_back(IntList{addedFaceCenterIdx, nodeIdxCurr, nodeIdxNext});
                            }
                            if (KMaxFace_.find(faceID) != KMaxFace_.end())
                            {
                                KMaxFaceNodes_.push_back(IntList{addedFaceCenterIdx, nodeIdxCurr, nodeIdxNext});
                            }
                        }
                    }
                    for (size_t j = 0; j < faceI.size(); j++)
                    {
                        const int nodeIdxCurr = faceI[j];
                        const int nodeIdxNext = faceI[j == faceI.size() - 1 ? 0 : j + 1];
                        decomposedCellNodes_.emplace_back(IntList{addedCellCenterIdx, addedFaceCenterIdx, nodeIdxCurr, nodeIdxNext});

                    }
                    continue;
                }
                else
                {
                    // not decomposed K face
                    IntList newCellNodes(faceI.size() + 1);
                    newCellNodes[faceI.size()] = addedCellCenterIdx;
                    for (size_t ni = 0; ni < faceI.size(); ni++)
                    {
                        newCellNodes[ni] = faceI[ni];
                    }
                    decomposedCellNodes_.emplace_back(newCellNodes);
                }
            }
            // add decomposed I J face cell
            IntList &cellKNodes = cellKNodes_[ci];
            decomposedCellNodes_.emplace_back(IntList{cellKNodes[0], cellKNodes[1], cellKNodes[5], cellKNodes[4], addedCellCenterIdx});
            decomposedCellNodes_.emplace_back(IntList{cellKNodes[1], cellKNodes[2], cellKNodes[6], cellKNodes[5], addedCellCenterIdx});
            decomposedCellNodes_.emplace_back(IntList{cellKNodes[2], cellKNodes[3], cellKNodes[7], cellKNodes[6], addedCellCenterIdx});
            decomposedCellNodes_.emplace_back(IntList{cellKNodes[3], cellKNodes[0], cellKNodes[4], cellKNodes[7], addedCellCenterIdx});
        }
        else // not decompose
        {
            const IntList &hexNodes = cellKNodes_[ci];
            decomposedCellNodes_.emplace_back(IntList{hexNodes[0], hexNodes[1], hexNodes[2], hexNodes[3], hexNodes[4], hexNodes[5], hexNodes[6], hexNodes[7]});

        }
    }
    // check nodeCell size;
    for (size_t i = 0; i < decomposedCellNodes_.size(); i++)
    {
        if (decomposedCellNodes_[i].size() != 4 && decomposedCellNodes_[i].size() != 5 && decomposedCellNodes_[i].size() != 8)
        {
            std::cerr << "Error: decomposedCellNodes_ size not match! "
                      << "Suport size= 4(tetrahedron) 5(pyramid) 8(hexaherdron)" << std::endl;
            std::abort();
        }
    }
    decompReservoirCellNum_ = int(decomposedCellNodes_.size());
    clock_t endT = std::clock();
    if (perfFlag_)
    {
        std::cout << "Cell decompose Time: " << double(endT - startT) / CLOCKS_PER_SEC << " s" << std::endl;
    }
}

void GRDECLGrid::genExtendCell(
    PointList &nodeCoord,
    IntListList &cellNodes,
    IntList &pointData,
    IntList &cellData)
{   
    // gen fault pointData
    pointData.resize(nodeCoord.size(), 0);
    {   
        for (Map<int, int>::iterator it = faceKValue_.begin(); it != faceKValue_.end(); it++)
        {
            const int faceIdx = it->first;
            const Face &faceI = faces_[faceIdx];
            for (size_t ni = 0; ni < faceI.size(); ni++)
            {
                const int nodeIdx = faceI[ni];
                pointData[nodeIdx] = it->second + 10;
            }
        }
        for (size_t fi = 0; fi < faultFaces_.size(); fi++)
        {
            const int faceIdx = faultFaces_[fi];
            const Face &faceI = faces_[faceIdx];
            for (size_t ni = 0; ni < faceI.size(); ni++)
            {
                const int nodeIdx = faceI[ni];
                pointData[nodeIdx] = 1;
            }
        }

    }
    // gen cellData
    cellData.resize(cellNodes.size(), 1);
    const int baseCellNum = int(cellNodes.size());
    // gen overburden cell
    IntSet IMaxNodeSet;
    IntSet IMinNodeSet;
    IntSet JMaxNodeSet;
    IntSet JMinNodeSet;
    IntSet KMaxNodeSet;
    IntSet KMinNodeSet;

    const size_t orgIMaxFaceSize = IMaxFaceNodes_.size();
    IntSet IMaxUnderFaceSet;
    IntSet IMaxOverFaceSet;
    const size_t orgIMinFaceSize = IMinFaceNodes_.size();
    IntSet IMinUnderFaceSet;
    IntSet IMinOverFaceSet;
    const size_t orgJMaxFaceSize = JMaxFaceNodes_.size();
    IntSet JMaxUnderFaceSet;
    IntSet JMaxOverFaceSet;
    const size_t orgJMinFaceSize = JMinFaceNodes_.size();
    IntSet JMinUnderFaceSet;
    IntSet JMinOverFaceSet;

    for(size_t i = 0; i < IMaxFaceNodes_.size(); i++)
    {
        for(size_t j = 0; j < IMaxFaceNodes_[i].size(); j++)
        {
            IMaxNodeSet.insert(IMaxFaceNodes_[i][j]);
        }
    }
    for(size_t i = 0; i < IMinFaceNodes_.size(); i++)
    {
        for(size_t j = 0; j < IMinFaceNodes_[i].size(); j++)
        {
            IMinNodeSet.insert(IMinFaceNodes_[i][j]);
        }
    }
    for(size_t i = 0; i < JMaxFaceNodes_.size(); i++)
    {
        for(size_t j = 0; j < JMaxFaceNodes_[i].size(); j++)
        {
            JMaxNodeSet.insert(JMaxFaceNodes_[i][j]);
        }
    }
    for(size_t i = 0; i < JMinFaceNodes_.size(); i++)
    {
        for(size_t j = 0; j < JMinFaceNodes_[i].size(); j++)
        {
            JMinNodeSet.insert(JMinFaceNodes_[i][j]);
        }
    }
    for(size_t i = 0; i < KMaxFaceNodes_.size(); i++)
    {
        for(size_t j = 0; j < KMaxFaceNodes_[i].size(); j++)
        {
            KMaxNodeSet.insert(KMaxFaceNodes_[i][j]);
        }
    }
    for(size_t i = 0; i < KMinFaceNodes_.size(); i++)
    {
        for(size_t j = 0; j < KMinFaceNodes_[i].size(); j++)
        {
            KMinNodeSet.insert(KMinFaceNodes_[i][j]);
        }
    }

    IntSet IMaxJMaxOverCornerNodes;
    IntSet IMaxJMaxUnderCornerNodes;
    IntSet IMaxJMinOverCornerNodes;
    IntSet IMaxJMinUnderCornerNodes;
    IntSet IMinJMaxOverCornerNodes;
    IntSet IMinJMaxUnderCornerNodes;
    IntSet IMinJMinOverCornerNodes;
    IntSet IMinJMinUnderCornerNodes;
    // extended K+ cell
    Map<int,int> nodeMap_;
    PointList basePoints;
    const double zPlus = std::abs(extendZMax_ - boundingBox_.second.z) / double(ZExtendNum_);
    int pCounter = 0;
    for(size_t fi = 0; fi < KMaxFaceNodes_.size(); fi++)
    {
        IntList& faceNodeIdx = KMaxFaceNodes_[fi];
        for(size_t ni = 0; ni < faceNodeIdx.size(); ni++)
        {
            if(nodeMap_.find(faceNodeIdx[ni]) == nodeMap_.end())
            {
                nodeMap_.emplace(faceNodeIdx[ni],pCounter);
                pCounter++;
                basePoints.push_back(nodeCoord[faceNodeIdx[ni]]);
            }
        }
    }
    int baseCoordSize = int(nodeCoord.size());
    for(int zi = 0 ; zi < ZExtendNum_; zi++)
    {   
        for(size_t pi = 0; pi < basePoints.size(); pi++)
        {
            nodeCoord.push_back(Point(basePoints[pi].x,basePoints[pi].y,basePoints[pi].z + zPlus * (zi + 1)));
            pointData.push_back(0);
        }
        // add overBurdenCellNodes_
        for(size_t ci = 0; ci < KMaxFaceNodes_.size(); ci++)
        {
            IntList& faceNodeIdx = KMaxFaceNodes_[ci];
            IntList newCellNodes = faceNodeIdx;
            IntList IMaxNodes;
            IntList IMinNodes;
            IntList JMaxNodes;
            IntList JMinNodes;
            for(size_t ni = 0; ni < faceNodeIdx.size(); ni++)
            {
                newCellNodes[ni] = (zi == 0) ? faceNodeIdx[ni] : baseCoordSize + (zi-1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[ni]];
                newCellNodes.push_back(int(baseCoordSize + zi * int(basePoints.size()) + nodeMap_[faceNodeIdx[ni]]));
                if(IMaxNodeSet.find(faceNodeIdx[ni]) != IMaxNodeSet.end())
                {
                    IMaxNodes.push_back(int(ni));
                }
                if (IMinNodeSet.find(faceNodeIdx[ni]) != IMinNodeSet.end())
                {
                    IMinNodes.push_back(int(ni));
                }
                if (JMaxNodeSet.find(faceNodeIdx[ni]) != JMaxNodeSet.end())
                {
                    JMaxNodes.push_back(int(ni));
                }
                if (JMinNodeSet.find(faceNodeIdx[ni]) != JMinNodeSet.end())
                {
                    JMinNodes.push_back(int(ni));
                }
                // for 4 corner nodes
                if (IMaxNodeSet.find(faceNodeIdx[ni]) != IMaxNodeSet.end() &&
                    JMaxNodeSet.find(faceNodeIdx[ni]) != JMaxNodeSet.end())
                {
                    IMaxJMaxOverCornerNodes.insert(newCellNodes[ni]);
                    IMaxJMaxOverCornerNodes.insert(newCellNodes[ni + faceNodeIdx.size()]);
                }
                if (IMinNodeSet.find(faceNodeIdx[ni]) != IMinNodeSet.end() &&
                    JMaxNodeSet.find(faceNodeIdx[ni]) != JMaxNodeSet.end())
                {
                    IMinJMaxOverCornerNodes.insert(newCellNodes[ni]);
                    IMinJMaxOverCornerNodes.insert(newCellNodes[ni + faceNodeIdx.size()]);
                }
                if (IMinNodeSet.find(faceNodeIdx[ni]) != IMinNodeSet.end() &&
                    JMinNodeSet.find(faceNodeIdx[ni]) != JMinNodeSet.end())
                {
                    IMinJMinOverCornerNodes.insert(newCellNodes[ni]);
                    IMinJMinOverCornerNodes.insert(newCellNodes[ni + faceNodeIdx.size()]);
                }
                if (IMaxNodeSet.find(faceNodeIdx[ni]) != IMaxNodeSet.end() &&
                    JMinNodeSet.find(faceNodeIdx[ni]) != JMinNodeSet.end())
                {
                    IMaxJMinOverCornerNodes.insert(newCellNodes[ni]);
                    IMaxJMinOverCornerNodes.insert(newCellNodes[ni + faceNodeIdx.size()]);
                }
            }
            cellNodes.push_back(newCellNodes);
            cellData.push_back(3);
            if (IMaxNodes.size() == 2)
            {
                if (IMaxNodes[0] > IMaxNodes[1])
                {
                    std::swap(IMaxNodes[0], IMaxNodes[1]);
                }
                IntList newXMaxFaceNodes;
                newXMaxFaceNodes.push_back((zi == 0) ? faceNodeIdx[IMaxNodes[0]] : baseCoordSize + (zi - 1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[IMaxNodes[0]]]);
                newXMaxFaceNodes.push_back(int(baseCoordSize + zi * int(basePoints.size()) + nodeMap_[faceNodeIdx[IMaxNodes[0]]]));
                newXMaxFaceNodes.push_back(int(baseCoordSize + zi * int(basePoints.size()) + nodeMap_[faceNodeIdx[IMaxNodes[1]]]));
                newXMaxFaceNodes.push_back((zi == 0) ? faceNodeIdx[IMaxNodes[1]] : baseCoordSize + (zi - 1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[IMaxNodes[1]]]);
                IMaxOverFaceSet.insert(int(IMaxFaceNodes_.size()));
                IMaxFaceNodes_.push_back(newXMaxFaceNodes);
            }
            if (IMinNodes.size() == 2)
            {
                if (IMinNodes[0] > IMinNodes[1])
                {
                    std::swap(IMinNodes[0], IMinNodes[1]);
                }
                IntList newXMinFaceNodes;
                newXMinFaceNodes.push_back((zi == 0) ? faceNodeIdx[IMinNodes[0]] : baseCoordSize + (zi - 1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[IMinNodes[0]]]);
                newXMinFaceNodes.push_back(int(baseCoordSize + zi * int(basePoints.size()) + nodeMap_[faceNodeIdx[IMinNodes[0]]]));
                newXMinFaceNodes.push_back(int(baseCoordSize + zi * int(basePoints.size()) + nodeMap_[faceNodeIdx[IMinNodes[1]]]));
                newXMinFaceNodes.push_back((zi == 0) ? faceNodeIdx[IMinNodes[1]] : baseCoordSize + (zi - 1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[IMinNodes[1]]]);
                IMinOverFaceSet.insert(int(IMinFaceNodes_.size()));
                IMinFaceNodes_.push_back(newXMinFaceNodes);
            }
            if (JMaxNodes.size() == 2)
            {
                if (JMaxNodes[0] > JMaxNodes[1])
                {
                    std::swap(JMaxNodes[0], JMaxNodes[1]);
                }
                IntList newYMaxFaceNodes;
                newYMaxFaceNodes.push_back((zi == 0) ? faceNodeIdx[JMaxNodes[0]] : baseCoordSize + (zi - 1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMaxNodes[0]]]);
                newYMaxFaceNodes.push_back(int(baseCoordSize + zi * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMaxNodes[0]]]));
                newYMaxFaceNodes.push_back(int(baseCoordSize + zi * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMaxNodes[1]]]));
                newYMaxFaceNodes.push_back((zi == 0) ? faceNodeIdx[JMaxNodes[1]] : baseCoordSize + (zi - 1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMaxNodes[1]]]);
                JMaxOverFaceSet.insert(int(JMaxFaceNodes_.size()));
                JMaxFaceNodes_.push_back(newYMaxFaceNodes);
            }
            if (JMinNodes.size() == 2)
            {
                if (JMinNodes[0] > JMinNodes[1])
                {
                    std::swap(JMinNodes[0], JMinNodes[1]);
                }
                IntList newYMinFaceNodes;
                newYMinFaceNodes.push_back((zi == 0) ? faceNodeIdx[JMinNodes[0]] : baseCoordSize + (zi - 1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMinNodes[0]]]);
                newYMinFaceNodes.push_back(int(baseCoordSize + zi * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMinNodes[0]]]));
                newYMinFaceNodes.push_back(int(baseCoordSize + zi * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMinNodes[1]]]));
                newYMinFaceNodes.push_back((zi == 0) ? faceNodeIdx[JMinNodes[1]] : baseCoordSize + (zi - 1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMinNodes[1]]]);
                JMinOverFaceSet.insert(int(JMinFaceNodes_.size()));
                JMinFaceNodes_.push_back(newYMinFaceNodes);
            }
        }
    }
    // extended K- cell
    //IntSet KMinNodeSet;
    nodeMap_.clear();
    basePoints.clear();
    const double zMinus = std::abs((boundingBox_.first.z - extendZMin_) / double(ZExtendNum_));
    pCounter = 0;
    for(size_t fi = 0; fi < KMinFaceNodes_.size(); fi++)
    {
        IntList& faceNodeIdx = KMinFaceNodes_[fi];
        for(size_t ni = 0; ni < faceNodeIdx.size(); ni++)
        {
            if(nodeMap_.find(faceNodeIdx[ni]) == nodeMap_.end())
            {
                nodeMap_.emplace(faceNodeIdx[ni],pCounter);
                pCounter++;
                basePoints.push_back(nodeCoord[faceNodeIdx[ni]]);
            }
        }
    }
    baseCoordSize = int(nodeCoord.size());
    for(int zi = 0 ; zi < ZExtendNum_; zi++)
    {   
        for(size_t pi = 0; pi < basePoints.size(); pi++)
        {
            nodeCoord.push_back(Point(basePoints[pi].x,basePoints[pi].y,basePoints[pi].z - zMinus * (zi + 1)));
            pointData.push_back(0);
        }
        // add underBurdenCellNodes_
        for(size_t ci = 0; ci < KMinFaceNodes_.size(); ci++)
        {
            IntList& faceNodeIdx = KMinFaceNodes_[ci];
            IntList newCellNodes = faceNodeIdx;
            IntList IMaxNodes;
            IntList IMinNodes;
            IntList JMaxNodes;
            IntList JMinNodes;
            for(size_t ni = 0; ni < faceNodeIdx.size(); ni++)
            {
                newCellNodes[ni] = (zi == 0) ? faceNodeIdx[ni] : baseCoordSize + (zi-1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[ni]];
                newCellNodes.push_back(int(baseCoordSize + zi * int(basePoints.size()) + nodeMap_[faceNodeIdx[ni]]));
                if(IMaxNodeSet.find(faceNodeIdx[ni]) != IMaxNodeSet.end())
                {
                    IMaxNodes.push_back(int(ni));
                }
                if (IMinNodeSet.find(faceNodeIdx[ni]) != IMinNodeSet.end())
                {
                    IMinNodes.push_back(int(ni));
                }
                if (JMaxNodeSet.find(faceNodeIdx[ni]) != JMaxNodeSet.end())
                {
                    JMaxNodes.push_back(int(ni));
                }
                if (JMinNodeSet.find(faceNodeIdx[ni]) != JMinNodeSet.end())
                {
                    JMinNodes.push_back(int(ni));
                }
                // for 4 corner nodes
                if (IMaxNodeSet.find(faceNodeIdx[ni]) != IMaxNodeSet.end() &&
                    JMaxNodeSet.find(faceNodeIdx[ni]) != JMaxNodeSet.end())
                {
                    IMaxJMaxUnderCornerNodes.insert(newCellNodes[ni]);
                    IMaxJMaxUnderCornerNodes.insert(newCellNodes[ni + faceNodeIdx.size()]);
                }
                if (IMinNodeSet.find(faceNodeIdx[ni]) != IMinNodeSet.end() &&
                    JMaxNodeSet.find(faceNodeIdx[ni]) != JMaxNodeSet.end())
                {
                    IMinJMaxUnderCornerNodes.insert(newCellNodes[ni]);
                    IMinJMaxUnderCornerNodes.insert(newCellNodes[ni + faceNodeIdx.size()]);
                }
                if (IMinNodeSet.find(faceNodeIdx[ni]) != IMinNodeSet.end() &&
                    JMinNodeSet.find(faceNodeIdx[ni]) != JMinNodeSet.end())
                {
                    IMinJMinUnderCornerNodes.insert(newCellNodes[ni]);
                    IMinJMinUnderCornerNodes.insert(newCellNodes[ni + faceNodeIdx.size()]);
                }
                if (IMaxNodeSet.find(faceNodeIdx[ni]) != IMaxNodeSet.end() &&
                    JMinNodeSet.find(faceNodeIdx[ni]) != JMinNodeSet.end())
                {
                    IMaxJMinUnderCornerNodes.insert(newCellNodes[ni]);
                    IMaxJMinUnderCornerNodes.insert(newCellNodes[ni + faceNodeIdx.size()]);
                }
            }
            cellNodes.push_back(newCellNodes);
            cellData.push_back(5);
            //underBurdenCellNodes_.push_back(newCellNodes);
            if (IMaxNodes.size() == 2)
            {
                if (IMaxNodes[0] > IMaxNodes[1])
                {
                    std::swap(IMaxNodes[0], IMaxNodes[1]);
                }
                IntList newXMaxFaceNodes;
                newXMaxFaceNodes.push_back((zi == 0) ? faceNodeIdx[IMaxNodes[0]] : baseCoordSize + (zi - 1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[IMaxNodes[0]]]);
                newXMaxFaceNodes.push_back(int(baseCoordSize + zi * int(basePoints.size()) + nodeMap_[faceNodeIdx[IMaxNodes[0]]]));
                newXMaxFaceNodes.push_back(int(baseCoordSize + zi * int(basePoints.size()) + nodeMap_[faceNodeIdx[IMaxNodes[1]]]));
                newXMaxFaceNodes.push_back((zi == 0) ? faceNodeIdx[IMaxNodes[1]] : baseCoordSize + (zi - 1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[IMaxNodes[1]]]);
                IMaxUnderFaceSet.insert(int(IMaxFaceNodes_.size()));
                IMaxFaceNodes_.push_back(newXMaxFaceNodes);
            }
            if (IMinNodes.size() == 2)
            {
                if (IMinNodes[0] > IMinNodes[1])
                {
                    std::swap(IMinNodes[0], IMinNodes[1]);
                }
                IntList newXMinFaceNodes;
                newXMinFaceNodes.push_back((zi == 0) ? faceNodeIdx[IMinNodes[0]] : baseCoordSize + (zi - 1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[IMinNodes[0]]]);
                newXMinFaceNodes.push_back(int(baseCoordSize + zi * int(basePoints.size()) + nodeMap_[faceNodeIdx[IMinNodes[0]]]));
                newXMinFaceNodes.push_back(int(baseCoordSize + zi * int(basePoints.size()) + nodeMap_[faceNodeIdx[IMinNodes[1]]]));
                newXMinFaceNodes.push_back((zi == 0) ? faceNodeIdx[IMinNodes[1]] : baseCoordSize + (zi - 1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[IMinNodes[1]]]);
                IMinUnderFaceSet.insert(int(IMinFaceNodes_.size()));
                IMinFaceNodes_.push_back(newXMinFaceNodes);
            }
            if (JMaxNodes.size() == 2)
            {
                if (JMaxNodes[0] > JMaxNodes[1])
                {
                    std::swap(JMaxNodes[0], JMaxNodes[1]);
                }
                IntList newYMaxFaceNodes;
                newYMaxFaceNodes.push_back((zi == 0) ? faceNodeIdx[JMaxNodes[0]] : baseCoordSize + (zi - 1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMaxNodes[0]]]);
                newYMaxFaceNodes.push_back(int(baseCoordSize + zi * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMaxNodes[0]]]));
                newYMaxFaceNodes.push_back(int(baseCoordSize + zi * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMaxNodes[1]]]));
                newYMaxFaceNodes.push_back((zi == 0) ? faceNodeIdx[JMaxNodes[1]] : baseCoordSize + (zi - 1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMaxNodes[1]]]);
                JMaxUnderFaceSet.insert(int(JMaxFaceNodes_.size()));
                JMaxFaceNodes_.push_back(newYMaxFaceNodes);
            }
            if (JMinNodes.size() == 2)
            {
                if (JMinNodes[0] > JMinNodes[1])
                {
                    std::swap(JMinNodes[0], JMinNodes[1]);
                }
                IntList newYMinFaceNodes;
                newYMinFaceNodes.push_back((zi == 0) ? faceNodeIdx[JMinNodes[0]] : baseCoordSize + (zi - 1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMinNodes[0]]]);
                newYMinFaceNodes.push_back(int(baseCoordSize + zi * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMinNodes[0]]]));
                newYMinFaceNodes.push_back(int(baseCoordSize + zi * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMinNodes[1]]]));
                newYMinFaceNodes.push_back((zi == 0) ? faceNodeIdx[JMinNodes[1]] : baseCoordSize + (zi - 1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMinNodes[1]]]);
                JMinUnderFaceSet.insert(int(JMinFaceNodes_.size()));
                JMinFaceNodes_.push_back(newYMinFaceNodes);
            }  
        }
    }
    // extended I+ cell
    nodeMap_.clear();
    basePoints.clear();
    const double extendXDis_ = std::abs(boundingBox_.second.x - boundingBox_.first.x) * extendXFactor_;
    const double xPlus = extendXDis_ / double(XExtendNum_);
    pCounter = 0;
    for(size_t fi = 0; fi < IMaxFaceNodes_.size(); fi++)
    {
        IntList& faceNodeIdx = IMaxFaceNodes_[fi];
        for(size_t ni = 0; ni < faceNodeIdx.size(); ni++)
        {
            if(nodeMap_.find(faceNodeIdx[ni]) == nodeMap_.end())
            {
                nodeMap_.emplace(faceNodeIdx[ni],pCounter);
                pCounter++;
                basePoints.push_back(nodeCoord[faceNodeIdx[ni]]);
            }
        }
    }
    baseCoordSize = int(nodeCoord.size());
    for(int xi = 0 ; xi < XExtendNum_; xi++)
    {   
        for(size_t pi = 0; pi < basePoints.size(); pi++)
        {
            nodeCoord.push_back(Point(basePoints[pi].x + xPlus * (xi + 1),basePoints[pi].y,basePoints[pi].z));
            pointData.push_back(0);
        }
        // add IMaxFaceCellNodes_
        for(size_t ci = 0; ci < IMaxFaceNodes_.size(); ci++)
        {
            IntList& faceNodeIdx = IMaxFaceNodes_[ci];
            IntList newCellNodes = faceNodeIdx;
            IntList JMaxNodes;
            IntList JMinNodes;
            IntList JMaxOverCornerNodes;
            IntList JMaxUnderCornerNodes;
            IntList JMinOverCornerNodes;
            IntList JMinUnderCornerNodes;
            for(size_t ni = 0; ni < faceNodeIdx.size(); ni++)
            {
                newCellNodes[ni] = (xi == 0) ? faceNodeIdx[ni] : baseCoordSize + (xi-1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[ni]];
                newCellNodes.push_back(int(baseCoordSize + xi * int(basePoints.size()) + nodeMap_[faceNodeIdx[ni]]));
                if (JMaxNodeSet.find(faceNodeIdx[ni]) != JMaxNodeSet.end())
                {
                    JMaxNodes.push_back(int(ni));
                }
                if (JMinNodeSet.find(faceNodeIdx[ni]) != JMinNodeSet.end())
                {
                    JMinNodes.push_back(int(ni));
                }
                if (IMaxJMaxOverCornerNodes.find(faceNodeIdx[ni]) != IMaxJMaxOverCornerNodes.end() ||
                    IMinJMaxOverCornerNodes.find(faceNodeIdx[ni]) != IMinJMaxOverCornerNodes.end())
                {
                    JMaxOverCornerNodes.push_back(int(ni));
                }
                if (IMinJMinOverCornerNodes.find(faceNodeIdx[ni]) != IMinJMinOverCornerNodes.end() ||
                    IMaxJMinOverCornerNodes.find(faceNodeIdx[ni]) != IMaxJMinOverCornerNodes.end())
                {
                    JMinOverCornerNodes.push_back(int(ni));
                }
                if (IMaxJMaxUnderCornerNodes.find(faceNodeIdx[ni]) != IMaxJMaxUnderCornerNodes.end() ||
                    IMinJMaxUnderCornerNodes.find(faceNodeIdx[ni]) != IMinJMaxUnderCornerNodes.end())
                {
                    JMaxUnderCornerNodes.push_back(int(ni));
                }
                if (IMinJMinUnderCornerNodes.find(faceNodeIdx[ni]) != IMinJMinUnderCornerNodes.end() ||
                    IMaxJMinUnderCornerNodes.find(faceNodeIdx[ni]) != IMaxJMinUnderCornerNodes.end())
                {
                    JMinUnderCornerNodes.push_back(int(ni));
                }
            }
            //IMaxFaceCellNodes_.push_back(newCellNodes);
            cellNodes.push_back(newCellNodes);
            if (ci < orgIMaxFaceSize) // sideburden
            {
                cellData.push_back(2);
            }
            else if (IMaxOverFaceSet.find(int(ci)) != IMaxOverFaceSet.end()) // over
            {
                cellData.push_back(4);
            }
            else if (IMaxUnderFaceSet.find(int(ci)) != IMaxUnderFaceSet.end()) // under
            {
                cellData.push_back(6);
            }
            else
            {
                cellData.push_back(2);
            }

            if (JMaxNodes.size() == 2)
            {
                if (JMaxNodes[0] > JMaxNodes[1])
                {
                    std::swap(JMaxNodes[0], JMaxNodes[1]);
                }
                IntList newJMaxFaceNodes;
                newJMaxFaceNodes.push_back((xi == 0) ? faceNodeIdx[JMaxNodes[0]] : baseCoordSize + (xi - 1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMaxNodes[0]]]);
                newJMaxFaceNodes.push_back(int(baseCoordSize + xi * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMaxNodes[0]]]));
                newJMaxFaceNodes.push_back(int(baseCoordSize + xi * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMaxNodes[1]]]));
                newJMaxFaceNodes.push_back((xi == 0) ? faceNodeIdx[JMaxNodes[1]] : baseCoordSize + (xi - 1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMaxNodes[1]]]);
                JMaxFaceNodes_.push_back(newJMaxFaceNodes);
            }
            if (JMinNodes.size() == 2)
            {
                if (JMinNodes[0] > JMinNodes[1])
                {
                    std::swap(JMinNodes[0], JMinNodes[1]);
                }
                IntList newJMinFaceNodes;
                newJMinFaceNodes.push_back((xi == 0) ? faceNodeIdx[JMinNodes[0]] : baseCoordSize + (xi - 1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMinNodes[0]]]);
                newJMinFaceNodes.push_back(int(baseCoordSize + xi * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMinNodes[0]]]));
                newJMinFaceNodes.push_back(int(baseCoordSize + xi * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMinNodes[1]]]));
                newJMinFaceNodes.push_back((xi == 0) ? faceNodeIdx[JMinNodes[1]] : baseCoordSize + (xi - 1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMinNodes[1]]]);
                JMinFaceNodes_.push_back(newJMinFaceNodes);
            }
            if (JMaxOverCornerNodes.size() == 2)
            {
                if (JMaxOverCornerNodes[0] > JMaxOverCornerNodes[1])
                {
                    std::swap(JMaxOverCornerNodes[0], JMaxOverCornerNodes[1]);
                }
                IntList newJMaxKMinFaceNodes;
                newJMaxKMinFaceNodes.push_back((xi == 0) ? faceNodeIdx[JMaxOverCornerNodes[0]] : baseCoordSize + (xi - 1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMaxOverCornerNodes[0]]]);
                newJMaxKMinFaceNodes.push_back(int(baseCoordSize + xi * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMaxOverCornerNodes[0]]]));
                newJMaxKMinFaceNodes.push_back(int(baseCoordSize + xi * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMaxOverCornerNodes[1]]]));
                newJMaxKMinFaceNodes.push_back((xi == 0) ? faceNodeIdx[JMaxOverCornerNodes[1]] : baseCoordSize + (xi - 1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMaxOverCornerNodes[1]]]);
                JMaxOverFaceSet.insert(int(JMaxFaceNodes_.size()));
                JMaxFaceNodes_.push_back(newJMaxKMinFaceNodes);
            }
            if (JMinOverCornerNodes.size() == 2)
            {
                if (JMinOverCornerNodes[0] > JMinOverCornerNodes[1])
                {
                    std::swap(JMinOverCornerNodes[0], JMinOverCornerNodes[1]);
                }
                IntList newJMinKMaxFaceNodes;
                newJMinKMaxFaceNodes.push_back((xi == 0) ? faceNodeIdx[JMinOverCornerNodes[0]] : baseCoordSize + (xi - 1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMinOverCornerNodes[0]]]);
                newJMinKMaxFaceNodes.push_back(int(baseCoordSize + xi * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMinOverCornerNodes[0]]]));
                newJMinKMaxFaceNodes.push_back(int(baseCoordSize + xi * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMinOverCornerNodes[1]]]));
                newJMinKMaxFaceNodes.push_back((xi == 0) ? faceNodeIdx[JMinOverCornerNodes[1]] : baseCoordSize + (xi - 1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMinOverCornerNodes[1]]]);
                JMinOverFaceSet.insert(int(JMinFaceNodes_.size()));
                JMinFaceNodes_.push_back(newJMinKMaxFaceNodes);
            }
            if (JMaxUnderCornerNodes.size() == 2)
            {
                if (JMaxUnderCornerNodes[0] > JMaxUnderCornerNodes[1])
                {
                    std::swap(JMaxUnderCornerNodes[0], JMaxUnderCornerNodes[1]);
                }
                IntList newJMaxKMaxFaceNodes;
                newJMaxKMaxFaceNodes.push_back((xi == 0) ? faceNodeIdx[JMaxUnderCornerNodes[0]] : baseCoordSize + (xi - 1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMaxUnderCornerNodes[0]]]);
                newJMaxKMaxFaceNodes.push_back(int(baseCoordSize + xi * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMaxUnderCornerNodes[0]]]));
                newJMaxKMaxFaceNodes.push_back(int(baseCoordSize + xi * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMaxUnderCornerNodes[1]]]));
                newJMaxKMaxFaceNodes.push_back((xi == 0) ? faceNodeIdx[JMaxUnderCornerNodes[1]] : baseCoordSize + (xi - 1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMaxUnderCornerNodes[1]]]);
                JMaxUnderFaceSet.insert(int(JMaxFaceNodes_.size()));
                JMaxFaceNodes_.push_back(newJMaxKMaxFaceNodes);
            }
            if (JMinUnderCornerNodes.size() == 2)
            {
                if (JMinUnderCornerNodes[0] > JMinUnderCornerNodes[1])
                {
                    std::swap(JMinUnderCornerNodes[0], JMinUnderCornerNodes[1]);
                }
                IntList newJMinKMinFaceNodes;
                newJMinKMinFaceNodes.push_back((xi == 0) ? faceNodeIdx[JMinUnderCornerNodes[0]] : baseCoordSize + (xi - 1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMinUnderCornerNodes[0]]]);
                newJMinKMinFaceNodes.push_back(int(baseCoordSize + xi * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMinUnderCornerNodes[0]]]));
                newJMinKMinFaceNodes.push_back(int(baseCoordSize + xi * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMinUnderCornerNodes[1]]]));
                newJMinKMinFaceNodes.push_back((xi == 0) ? faceNodeIdx[JMinUnderCornerNodes[1]] : baseCoordSize + (xi - 1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMinUnderCornerNodes[1]]]);
                JMinUnderFaceSet.insert(int(JMinFaceNodes_.size()));
                JMinFaceNodes_.push_back(newJMinKMinFaceNodes);
            }
        }
    }
    // extended I- cell
    nodeMap_.clear();
    basePoints.clear();
    const double xMinus = extendXDis_ / double(XExtendNum_);
    pCounter = 0;
    for(size_t fi = 0; fi < IMinFaceNodes_.size(); fi++)
    {
        IntList& faceNodeIdx = IMinFaceNodes_[fi];
        for(size_t ni = 0; ni < faceNodeIdx.size(); ni++)
        {
            if(nodeMap_.find(faceNodeIdx[ni]) == nodeMap_.end())
            {
                nodeMap_.emplace(faceNodeIdx[ni],pCounter);
                pCounter++;
                basePoints.push_back(nodeCoord[faceNodeIdx[ni]]);
            }
        }
    }
    baseCoordSize = int(nodeCoord.size());
    for(int xi = 0 ; xi < XExtendNum_; xi++)
    {   
        for(size_t pi = 0; pi < basePoints.size(); pi++)
        {
            nodeCoord.push_back(Point(basePoints[pi].x - xMinus * (xi + 1),basePoints[pi].y,basePoints[pi].z));
            pointData.push_back(0);
        }
        // add IMinFaceCellNodes_
        for(size_t ci = 0; ci < IMinFaceNodes_.size(); ci++)
        {
            IntList& faceNodeIdx = IMinFaceNodes_[ci];
            IntList newCellNodes = faceNodeIdx;
            IntList JMaxNodes;
            IntList JMinNodes;
            IntList JMaxOverCornerNodes;
            IntList JMaxUnderCornerNodes;
            IntList JMinOverCornerNodes;
            IntList JMinUnderCornerNodes;
            for(size_t ni = 0; ni < faceNodeIdx.size(); ni++)
            {
                newCellNodes[ni] = (xi == 0) ? faceNodeIdx[ni] : baseCoordSize + (xi-1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[ni]];
                newCellNodes.push_back(int(baseCoordSize + xi * int(basePoints.size()) + nodeMap_[faceNodeIdx[ni]]));
                if (JMaxNodeSet.find(faceNodeIdx[ni]) != JMaxNodeSet.end())
                {
                    JMaxNodes.push_back(int(ni));
                }
                if (JMinNodeSet.find(faceNodeIdx[ni]) != JMinNodeSet.end())
                {
                    JMinNodes.push_back(int(ni));
                }
                if (IMaxJMaxOverCornerNodes.find(faceNodeIdx[ni]) != IMaxJMaxOverCornerNodes.end() ||
                    IMinJMaxOverCornerNodes.find(faceNodeIdx[ni]) != IMinJMaxOverCornerNodes.end())
                {
                    JMaxOverCornerNodes.push_back(int(ni));
                }
                if (IMinJMinOverCornerNodes.find(faceNodeIdx[ni]) != IMinJMinOverCornerNodes.end() ||
                    IMaxJMinOverCornerNodes.find(faceNodeIdx[ni]) != IMaxJMinOverCornerNodes.end())
                {
                    JMinOverCornerNodes.push_back(int(ni));
                }
                if (IMaxJMaxUnderCornerNodes.find(faceNodeIdx[ni]) != IMaxJMaxUnderCornerNodes.end() ||
                    IMinJMaxUnderCornerNodes.find(faceNodeIdx[ni]) != IMinJMaxUnderCornerNodes.end())
                {
                    JMaxUnderCornerNodes.push_back(int(ni));
                }
                if (IMinJMinUnderCornerNodes.find(faceNodeIdx[ni]) != IMinJMinUnderCornerNodes.end() ||
                    IMaxJMinUnderCornerNodes.find(faceNodeIdx[ni]) != IMaxJMinUnderCornerNodes.end())
                {
                    JMinUnderCornerNodes.push_back(int(ni));
                }
            }
            cellNodes.push_back(newCellNodes);
            if (ci < orgIMinFaceSize) // sideburden
            {
                cellData.push_back(2);
            }
            else if (IMinOverFaceSet.find(int(ci)) != IMinOverFaceSet.end()) // over
            {
                cellData.push_back(4);
            }
            else if (IMinUnderFaceSet.find(int(ci)) != IMinUnderFaceSet.end()) // under
            {
                cellData.push_back(6);
            }
            else
            {
                cellData.push_back(2);
            }

            //IMinFaceCellNodes_.push_back(newCellNodes);
            if (JMaxNodes.size() == 2)
            {
                if (JMaxNodes[0] > JMaxNodes[1])
                {
                    std::swap(JMaxNodes[0], JMaxNodes[1]);
                }
                IntList newJMaxFaceNodes;
                newJMaxFaceNodes.push_back((xi == 0) ? faceNodeIdx[JMaxNodes[0]] : baseCoordSize + (xi - 1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMaxNodes[0]]]);
                newJMaxFaceNodes.push_back(int(baseCoordSize + xi * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMaxNodes[0]]]));
                newJMaxFaceNodes.push_back(int(baseCoordSize + xi * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMaxNodes[1]]]));
                newJMaxFaceNodes.push_back((xi == 0) ? faceNodeIdx[JMaxNodes[1]] : baseCoordSize + (xi - 1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMaxNodes[1]]]);
                JMaxFaceNodes_.push_back(newJMaxFaceNodes);
            }
            if (JMinNodes.size() == 2)
            {
                if (JMinNodes[0] > JMinNodes[1])
                {
                    std::swap(JMinNodes[0], JMinNodes[1]);
                }
                IntList newJMinFaceNodes;
                newJMinFaceNodes.push_back((xi == 0) ? faceNodeIdx[JMinNodes[0]] : baseCoordSize + (xi - 1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMinNodes[0]]]);
                newJMinFaceNodes.push_back(int(baseCoordSize + xi * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMinNodes[0]]]));
                newJMinFaceNodes.push_back(int(baseCoordSize + xi * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMinNodes[1]]]));
                newJMinFaceNodes.push_back((xi == 0) ? faceNodeIdx[JMinNodes[1]] : baseCoordSize + (xi - 1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMinNodes[1]]]);
                JMinFaceNodes_.push_back(newJMinFaceNodes);
            }
            if (JMaxOverCornerNodes.size() == 2)
            {
                if (JMaxOverCornerNodes[0] > JMaxOverCornerNodes[1])
                {
                    std::swap(JMaxOverCornerNodes[0], JMaxOverCornerNodes[1]);
                }
                IntList newJMaxKMinFaceNodes;
                newJMaxKMinFaceNodes.push_back((xi == 0) ? faceNodeIdx[JMaxOverCornerNodes[0]] : baseCoordSize + (xi - 1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMaxOverCornerNodes[0]]]);
                newJMaxKMinFaceNodes.push_back(int(baseCoordSize + xi * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMaxOverCornerNodes[0]]]));
                newJMaxKMinFaceNodes.push_back(int(baseCoordSize + xi * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMaxOverCornerNodes[1]]]));
                newJMaxKMinFaceNodes.push_back((xi == 0) ? faceNodeIdx[JMaxOverCornerNodes[1]] : baseCoordSize + (xi - 1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMaxOverCornerNodes[1]]]);
                JMaxOverFaceSet.insert(int(JMaxFaceNodes_.size()));
                JMaxFaceNodes_.push_back(newJMaxKMinFaceNodes);
            }
            if (JMinOverCornerNodes.size() == 2)
            {
                if (JMinOverCornerNodes[0] > JMinOverCornerNodes[1])
                {
                    std::swap(JMinOverCornerNodes[0], JMinOverCornerNodes[1]);
                }
                IntList newJMinKMaxFaceNodes;
                newJMinKMaxFaceNodes.push_back((xi == 0) ? faceNodeIdx[JMinOverCornerNodes[0]] : baseCoordSize + (xi - 1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMinOverCornerNodes[0]]]);
                newJMinKMaxFaceNodes.push_back(int(baseCoordSize + xi * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMinOverCornerNodes[0]]]));
                newJMinKMaxFaceNodes.push_back(int(baseCoordSize + xi * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMinOverCornerNodes[1]]]));
                newJMinKMaxFaceNodes.push_back((xi == 0) ? faceNodeIdx[JMinOverCornerNodes[1]] : baseCoordSize + (xi - 1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMinOverCornerNodes[1]]]);
                JMinOverFaceSet.insert(int(JMinFaceNodes_.size()));
                JMinFaceNodes_.push_back(newJMinKMaxFaceNodes);
            }
            if (JMaxUnderCornerNodes.size() == 2)
            {
                if (JMaxUnderCornerNodes[0] > JMaxUnderCornerNodes[1])
                {
                    std::swap(JMaxUnderCornerNodes[0], JMaxUnderCornerNodes[1]);
                }
                IntList newJMaxKMinFaceNodes;
                newJMaxKMinFaceNodes.push_back((xi == 0) ? faceNodeIdx[JMaxUnderCornerNodes[0]] : baseCoordSize + (xi - 1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMaxUnderCornerNodes[0]]]);
                newJMaxKMinFaceNodes.push_back(int(baseCoordSize + xi * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMaxUnderCornerNodes[0]]]));
                newJMaxKMinFaceNodes.push_back(int(baseCoordSize + xi * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMaxUnderCornerNodes[1]]]));
                newJMaxKMinFaceNodes.push_back((xi == 0) ? faceNodeIdx[JMaxUnderCornerNodes[1]] : baseCoordSize + (xi - 1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMaxUnderCornerNodes[1]]]);
                JMaxUnderFaceSet.insert(int(JMaxFaceNodes_.size()));
                JMaxFaceNodes_.push_back(newJMaxKMinFaceNodes);
            }
            if (JMinUnderCornerNodes.size() == 2)
            {
                if (JMinUnderCornerNodes[0] > JMinUnderCornerNodes[1])
                {
                    std::swap(JMinUnderCornerNodes[0], JMinUnderCornerNodes[1]);
                }
                IntList newJMinKMaxFaceNodes;
                newJMinKMaxFaceNodes.push_back((xi == 0) ? faceNodeIdx[JMinUnderCornerNodes[0]] : baseCoordSize + (xi - 1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMinUnderCornerNodes[0]]]);
                newJMinKMaxFaceNodes.push_back(int(baseCoordSize + xi * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMinUnderCornerNodes[0]]]));
                newJMinKMaxFaceNodes.push_back(int(baseCoordSize + xi * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMinUnderCornerNodes[1]]]));
                newJMinKMaxFaceNodes.push_back((xi == 0) ? faceNodeIdx[JMinUnderCornerNodes[1]] : baseCoordSize + (xi - 1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[JMinUnderCornerNodes[1]]]);
                JMinUnderFaceSet.insert(int(JMinFaceNodes_.size()));
                JMinFaceNodes_.push_back(newJMinKMaxFaceNodes);
            }
        }
    }

    // extended J+ cell
    nodeMap_.clear();
    basePoints.clear();
    const double extendYDis_ = std::abs(boundingBox_.second.y - boundingBox_.first.y) * extendYFactor_;
    const double yPlus = extendYDis_ / double(YExtendNum_);
    pCounter = 0;
    for(size_t fi = 0; fi < JMaxFaceNodes_.size(); fi++)
    {
        IntList& faceNodeIdx = JMaxFaceNodes_[fi];
        for(size_t ni = 0; ni < faceNodeIdx.size(); ni++)
        {
            if(nodeMap_.find(faceNodeIdx[ni]) == nodeMap_.end())
            {
                nodeMap_.emplace(faceNodeIdx[ni],pCounter);
                pCounter++;
                basePoints.push_back(nodeCoord[faceNodeIdx[ni]]);
            }
        }
    }
    baseCoordSize = int(nodeCoord.size());
    for(int yi = 0 ; yi < YExtendNum_; yi++)
    {   
        for(size_t pi = 0; pi < basePoints.size(); pi++)
        {
            nodeCoord.push_back(Point(basePoints[pi].x,basePoints[pi].y - yPlus * (yi + 1),basePoints[pi].z));
            pointData.push_back(0);
        }
        // add JMaxFaceCellNodes_
        for(size_t ci = 0; ci < JMaxFaceNodes_.size(); ci++)
        {
            IntList& faceNodeIdx = JMaxFaceNodes_[ci];
            IntList newCellNodes = faceNodeIdx;
            for(size_t ni = 0; ni < faceNodeIdx.size(); ni++)
            {
                newCellNodes[ni] = (yi == 0) ? faceNodeIdx[ni] : baseCoordSize + (yi-1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[ni]];
                newCellNodes.push_back(int(baseCoordSize + yi * int(basePoints.size()) + nodeMap_[faceNodeIdx[ni]]));
            }
            if (ci < orgJMaxFaceSize) // sideburden
            {
                cellData.push_back(2);
            }
            else if (JMaxOverFaceSet.find(int(ci)) != JMaxOverFaceSet.end()) // over
            {
                cellData.push_back(4);
            }
            else if (JMaxUnderFaceSet.find(int(ci)) != JMaxUnderFaceSet.end()) // under
            {
                cellData.push_back(6);
            }
            else
            {
                cellData.push_back(2);
            }
            cellNodes.push_back(newCellNodes);
            //JMaxFaceCellNodes_.push_back(newCellNodes);
        }
    }
    //extended J- cell
    nodeMap_.clear();
    basePoints.clear();
    const double yMinus = extendYDis_ / double(YExtendNum_);
    pCounter = 0;
    for(size_t fi = 0; fi < JMinFaceNodes_.size(); fi++)
    {
        IntList& faceNodeIdx = JMinFaceNodes_[fi];
        for(size_t ni = 0; ni < faceNodeIdx.size(); ni++)
        {
            if(nodeMap_.find(faceNodeIdx[ni]) == nodeMap_.end())
            {
                nodeMap_.emplace(faceNodeIdx[ni],pCounter);
                pCounter++;
                basePoints.push_back(nodeCoord[faceNodeIdx[ni]]);
            }
        }
    }
    baseCoordSize = int(nodeCoord.size());
    for(int yi = 0 ; yi < YExtendNum_; yi++)
    {   
        for(size_t pi = 0; pi < basePoints.size(); pi++)
        {
            nodeCoord.push_back(Point(basePoints[pi].x,basePoints[pi].y + yMinus * (yi + 1),basePoints[pi].z));
            pointData.push_back(0);
        }
        // add JMinFaceCellNodes_
        for(size_t ci = 0; ci < JMinFaceNodes_.size(); ci++)
        {
            IntList& faceNodeIdx = JMinFaceNodes_[ci];
            IntList newCellNodes = faceNodeIdx;
            for(size_t ni = 0; ni < faceNodeIdx.size(); ni++)
            {
                newCellNodes[ni] = (yi == 0) ? faceNodeIdx[ni] : baseCoordSize + (yi-1) * int(basePoints.size()) + nodeMap_[faceNodeIdx[ni]];
                newCellNodes.push_back(int(baseCoordSize + yi * int(basePoints.size()) + nodeMap_[faceNodeIdx[ni]]));
            }
            cellNodes.push_back(newCellNodes);
            if (ci < orgJMinFaceSize) // sideburden
            {
                cellData.push_back(2);
            }
            else if (JMinOverFaceSet.find(int(ci)) != JMinOverFaceSet.end()) // over
            {
                cellData.push_back(4);
            }
            else if (JMinUnderFaceSet.find(int(ci)) != JMinUnderFaceSet.end()) // under
            {
                cellData.push_back(6);
            }
            else
            {
                cellData.push_back(2);
            }
            //JMinFaceCellNodes_.push_back(newCellNodes);
        }
    }
}

void GRDECLGrid::vtkCellDataWriter(const std::string &outputPath, const DoubleList &data, const DoubleList &data2) const
{
    if (data.size() != size_t(nCell_))
    {
        std::cerr << "Data size not match with points size!" << std::endl;
        std::abort();
    }

    std::ofstream outputFile(outputPath);

    if (!outputFile)
    {
        std::cerr << " Can not open: " << outputPath << std::endl;
        std::abort();
    }

    /**********Write header***********/
    outputFile << "# vtk DataFile Version 2.0\n";
    outputFile << "GRDECL Mesh\n";
    outputFile << "ASCII\n";
    outputFile << "DATASET UNSTRUCTURED_GRID\n";
    outputFile << "POINTS " << nodeCoord_.size() << " double"
               << "\n";
    /**********Write Points Data******************/
    outputFile << std::setprecision(11);
    //  "Writing fracture mesh point" << endl;
    for (size_t i = 0; i < nodeCoord_.size(); i++)
    {
        outputFile << nodeCoord_[i].x << "  ";
        outputFile << nodeCoord_[i].y << "  ";
        outputFile << nodeCoord_[i].z << "\n";
    }
    outputFile << "\n";
    /**********Write CEll************/
    // CELL header
    //  "Writing fracture mesh face" << endl;
    outputFile << "CELLS " << nCell_ << " " << nCell_ * 9 << "\n";
    for (size_t i = 0; i <cellKNodes_.size(); i++)
    {
        outputFile << "8 ";
        for (size_t j = 0; j <cellKNodes_[i].size(); j++)
        {
            outputFile <<cellKNodes_[i][j] << " ";
        }
        outputFile << "\n";
    }
    outputFile << "\n";

    // Writing cell types;
    outputFile << "CELL_TYPES " << nCell_ << '\n';
    for (int i = 0; i < nCell_; i++)
    {
        outputFile << "12\n";
    }
    outputFile << '\n';
    /**********Write CEll Data************/
    // Writing cell data;
    outputFile << "CELL_DATA " << nCell_ << '\n';
    outputFile << "SCALARS Data1YoungMode double 1\n";
    outputFile << "LOOKUP_TABLE Data1YoungMode\n";
    for (size_t i = 0; i < data.size(); i++)
    {
        outputFile << data[i] << '\n';
    }
    outputFile << '\n';

    outputFile << "SCALARS Data2Possion double 1\n";
    outputFile << "LOOKUP_TABLE Data2Possion\n";
    for (size_t i = 0; i < data2.size(); i++)
    {
        outputFile << data2[i] << '\n';
    }

}

void GRDECLGrid::vtkHexCellWriter(const std::string &outputPath)
{

    std::ofstream outputFile(outputPath);

    if (!outputFile)
    {
        std::cerr << " Can not open: " << outputPath << std::endl;
        std::abort();
    }

    /**********Write header***********/
    outputFile << "# vtk DataFile Version 2.0\n";
    outputFile << "GRDECL Mesh\n";
    outputFile << "ASCII\n";
    outputFile << "DATASET UNSTRUCTURED_GRID\n";
    outputFile << "POINTS " << allHexCoord_.size() << " double"
               << "\n";
    /**********Write Points Data******************/
    outputFile << std::setprecision(11);
    for (size_t i = 0; i < allHexCoord_.size(); i++)
    {
        outputFile << allHexCoord_[i].x << "  ";
        outputFile << allHexCoord_[i].y << "  ";
        outputFile << allHexCoord_[i].z << "\n";
    }
    outputFile << "\n";
    /**********Write CEll Data************/
    // CELL header
    outputFile << "CELLS " << allHexCellNodes_.size() << " " << allHexCellNodes_.size() * 9 << "\n";
    for (size_t i = 0; i < allHexCellNodes_.size(); i++)
    {
        outputFile << "8 ";
        for (size_t j = 0; j < allHexCellNodes_[i].size(); j++)
        {
            outputFile << allHexCellNodes_[i][j] << " ";
        }
        outputFile << "\n";
    }
    outputFile << "\n";

    // Writing cell types;
    outputFile << "CELL_TYPES " << allHexCellNodes_.size() << '\n';
    for (size_t i = 0; i < allHexCellNodes_.size(); i++)
    {
        outputFile << "12\n";
    }

    /**********Write POINT Data************/
    //prepare data 
    IntSet faultPoints; // Also: global to loacl Map
    for (size_t fi = 0; fi < faultFaces_.size(); fi++)
    {
        const int faceIdx = faultFaces_[fi];
        const Face &faceI = faces_[faceIdx];
        for (size_t ni = 0; ni < faceI.size(); ni++)
        {
            const int nodeIdx = faceI[ni];
            if (faultPoints.find(nodeIdx) == faultPoints.end())
            {
                faultPoints.insert(nodeIdx);
            }
        }
    }
    outputFile << "POINT_DATA " << allHexCoord_.size() << '\n';
    outputFile << "SCALARS FaultTag double 1\n";
    outputFile << "LOOKUP_TABLE FaultTag\n";
    for (size_t i = 0; i < allHexCoord_.size(); i++)
    {
        if (faultPoints.find(int(i)) != faultPoints.end())
        {
            outputFile << "1\n";
        }
        else
        {
            outputFile << "0\n";
        }
    }

    outputFile << '\n';

    // free data for memory test
    // allHexCoord_.resize(0);
    // allHexCellNodes_.resize(0);

}

void GRDECLGrid::vtkTetCellWriter(const std::string &outputPath)
{
    std::ofstream outputFile(outputPath);

    if (!outputFile)
    {
        std::cerr << " Can not open: " << outputPath << std::endl;
        std::abort();
    }

    /**********Write header***********/
    outputFile << "# vtk DataFile Version 2.0\n";
    outputFile << "GRDECL Mesh\n";
    outputFile << "ASCII\n";
    outputFile << "DATASET UNSTRUCTURED_GRID\n";
    outputFile << "POINTS " << allTetCoord_.size() << " double"
               << "\n";
    /**********Write Points Data******************/
    outputFile << std::setprecision(11);
    for (size_t i = 0; i < allTetCoord_.size(); i++)
    {
        outputFile << allTetCoord_[i].x << "  ";
        outputFile << allTetCoord_[i].y << "  ";
        outputFile << allTetCoord_[i].z << "\n";
    }
    outputFile << "\n";
    /**********Write CEll Data************/
    // CELL header
    //  "Writing fracture mesh face" << endl;
    outputFile << "CELLS " << allTetCellNodes_.size() << " " << allTetCellNodes_.size() * 5 << "\n";
    for (size_t i = 0; i < allTetCellNodes_.size(); i++)
    {
        outputFile << "4 ";
        for (size_t j = 0; j < allTetCellNodes_[i].size(); j++)
        {
            outputFile << allTetCellNodes_[i][j] << " ";
        }
        outputFile << "\n";
    }
    outputFile << "\n";

    // Writing cell types;
    outputFile << "CELL_TYPES " << allTetCellNodes_.size() << '\n';
    for (size_t i = 0; i < allTetCellNodes_.size(); i++)
    {
        outputFile << "10\n";
    }

    /**********Write POINT Data************/
    //prepare data 
    IntSet faultPoints; // Also: global to loacl Map
    for (size_t fi = 0; fi < faultFaces_.size(); fi++)
    {
        const int faceIdx = faultFaces_[fi];
        const Face &faceI = faces_[faceIdx];
        for (size_t ni = 0; ni < faceI.size(); ni++)
        {
            const int nodeIdx = faceI[ni];
            if (faultPoints.find(nodeIdx) == faultPoints.end())
            {
                faultPoints.insert(nodeIdx);
            }
        }
    }
    outputFile << "POINT_DATA " << allTetCoord_.size() << '\n';
    outputFile << "SCALARS FaultTag double 1\n";
    outputFile << "LOOKUP_TABLE FaultTag\n";
    for (size_t i = 0; i < allTetCoord_.size(); i++)
    {
        if (faultPoints.find(int(i)) != faultPoints.end())
        {
            outputFile << "1\n";
        }
        else
        {
            outputFile << "0\n";
        }
    }
    outputFile << '\n';

    // free data for memory test
    allTetCoord_.resize(0);
    allTetCellNodes_.resize(0);
}

void GRDECLGrid::vtkDecomposeCellWriter(const std::string &outputPath) const
{
    std::ofstream outputFile(outputPath);

    if (!outputFile)
    {
        std::cerr << " Can not open: " << outputPath << std::endl;
        std::abort();
    }

    /**********Write header***********/
    outputFile << "# vtk DataFile Version 2.0\n";
    outputFile << "GRDECL Mesh\n";
    outputFile << "ASCII\n";
    outputFile << "DATASET UNSTRUCTURED_GRID\n";
    outputFile << "POINTS " << decomposedCoord_.size() << " double"
               << "\n";
    /**********Write Points Data******************/
    outputFile << std::setprecision(11);
    for (size_t i = 0; i < decomposedCoord_.size(); i++)
    {
        outputFile << decomposedCoord_[i].x << "  ";
        outputFile << decomposedCoord_[i].y << "  ";
        outputFile << decomposedCoord_[i].z << "\n";
    }
    outputFile << "\n";
    /**********Write CEll Data************/
    // CELL header
    // loop cell for find sum
    size_t cellNodesCounter = 0;
    for (size_t i = 0; i < decomposedCellNodes_.size(); i++)
    {
        cellNodesCounter += decomposedCellNodes_[i].size();
    }
    int sizeTypeMap[8] = {-1, -1, -1, 10, 14, -1, -1, 12};
    // Writing cell data;
    outputFile << "CELLS " << decomposedCellNodes_.size() << " " << decomposedCellNodes_.size() + cellNodesCounter << "\n";
    for (size_t i = 0; i < decomposedCellNodes_.size(); i++)
    {
        outputFile << decomposedCellNodes_[i].size() << " ";
        for (size_t j = 0; j < decomposedCellNodes_[i].size(); j++)
        {
            outputFile << decomposedCellNodes_[i][j] << " ";
        }
        outputFile << "\n";
    }
    outputFile << "\n";

    // Writing cell types;
    outputFile << "CELL_TYPES " << decomposedCellNodes_.size() << '\n';
    for (size_t i = 0; i < decomposedCellNodes_.size(); i++)
    {
        outputFile << sizeTypeMap[decomposedCellNodes_[i].size()-1] << '\n';
    }

    /**********Write POINT Data************/
    //prepare data 
    IntSet faultPoints; // Also: global to loacl Map
    
    for (size_t fi = 0; fi < faultFaces_.size(); fi++)
    {
        const int faceIdx = faultFaces_[fi];
        const Face &faceI = faces_[faceIdx];
        for (size_t ni = 0; ni < faceI.size(); ni++)
        {
            const int nodeIdx = faceI[ni];
            if (faultPoints.find(nodeIdx) == faultPoints.end())
            {
                faultPoints.insert(nodeIdx);
            }
        }
    }
    outputFile << "POINT_DATA " << decomposedCoord_.size() << '\n';
    outputFile << "SCALARS FaultTag double 1\n";
    outputFile << "LOOKUP_TABLE FaultTag\n";
    for (size_t i = 0; i < decomposedCoord_.size(); i++)
    {
        if (faultPoints.find(int(i)) != faultPoints.end())
        {
            outputFile << "1\n";
        }
        else
        {
            outputFile << "0\n";
        }
    }
    outputFile << '\n';
}

void GRDECLGrid::vtkCellWriter(const std::string &outputPath) const
{
    std::ofstream outputFile(outputPath);

    if (!outputFile)
    {
        std::cerr << " Can not open: " << outputPath << std::endl;
        std::abort();
    }

    /**********Write header***********/
    outputFile << "# vtk DataFile Version 2.0\n";
    outputFile << "GRDECL Mesh\n";
    outputFile << "ASCII\n";
    outputFile << "DATASET UNSTRUCTURED_GRID\n";
    outputFile << "POINTS " << nodeCoord_.size() << " double"
               << "\n";
    /**********Write Points Data******************/
    //  "Writing mesh point" << endl;
    outputFile << std::setprecision(11);
    for (size_t i = 0; i < nodeCoord_.size(); i++)
    {
        outputFile << nodeCoord_[i].x << "  ";
        outputFile << nodeCoord_[i].y << "  ";
        outputFile << nodeCoord_[i].z << "\n";
    }
    outputFile << "\n";
    /**********Write CEll Data************/
    int counter = 0;
    for (size_t i = 0; i < faces_.size(); i++)
    {
        counter += int(faces_[i].size());
    }
    // CELL header
    //  "Writing mesh face" << endl;
    outputFile << "CELLS " << faces_.size() << " " << faces_.size() + counter << "\n";
    for (size_t i = 0; i < faces_.size(); i++)
    {
        outputFile <<faces_[i].size() <<" ";
        for (size_t j = 0; j < faces_[i].size(); j++)
        {
            outputFile << faces_[i][j] << " ";
        }
        outputFile << "\n";
    }
    outputFile << "\n";

    // Writing cell types;
    outputFile << "CELL_TYPES " << faces_.size() << '\n';
    for (size_t i = 0; i < faces_.size(); i++)
    {
        outputFile << "7\n";
    }
    outputFile << '\n';
}

void GRDECLGrid::vtkFaultFaceWriter(const std::string &outputPath) const
{
    std::ofstream outputFile(outputPath);

    if (!outputFile)
    {
        std::cerr << " Can not open: " << outputPath << std::endl;
        std::abort();
    }

    /**********Prepare Data***********/
    Map<int, int> addedPointsMap; // Also: global to loacl Map
    IntList localPointToGlobalMap;
    for (size_t fi = 0; fi < faultFaces_.size(); fi++)
    {
        const int faceIdx = faultFaces_[fi];
        const Face &faceI = faces_[faceIdx];
        for (size_t ni = 0; ni < faceI.size(); ni++)
        {
            const int nodeIdx = faceI[ni];
            if (addedPointsMap.find(nodeIdx) == addedPointsMap.end())
            {
                addedPointsMap.emplace(nodeIdx, int(localPointToGlobalMap.size()));
                localPointToGlobalMap.push_back(nodeIdx);
            }
        }
    }

    /**********Write header***********/
    outputFile << "# vtk DataFile Version 2.0\n";
    outputFile << "GRDECL Mesh\n";
    outputFile << "ASCII\n";
    outputFile << "DATASET UNSTRUCTURED_GRID\n";
    outputFile << "POINTS " << localPointToGlobalMap.size() << " double"
               << "\n";
    /**********Write Points Data******************/
    outputFile << std::setprecision(11);
    for (size_t i = 0; i < localPointToGlobalMap.size(); i++)
    {
        const int gIdx = localPointToGlobalMap[i];
        outputFile << nodeCoord_[gIdx].x << "  ";
        outputFile << nodeCoord_[gIdx].y << "  ";
        outputFile << nodeCoord_[gIdx].z << "\n";
    }
    outputFile << "\n";
    int counter = 0;
    for (size_t i = 0; i < faultFaces_.size(); i++)
    {
        counter += int(faces_[faultFaces_[i]].size());
    }

    outputFile << "CELLS " << faultFaces_.size() << " " << faultFaces_.size() + counter << "\n";
    for (size_t i = 0; i < faultFaces_.size(); i++)
    {
        const int faceIdx = faultFaces_[i];
        outputFile << faces_[faceIdx].size() << " ";
        for (size_t j = 0; j < faces_[faceIdx].size(); j++)
        {
            outputFile << addedPointsMap[faces_[faceIdx][j]] << " ";
        }
        outputFile << "\n";
    }
    outputFile << "\n";

    outputFile << "CELL_TYPES " << faultFaces_.size() << '\n';
    for (size_t i = 0; i < faultFaces_.size(); i++)
    {
        outputFile << "7\n";
    }
    outputFile << '\n';
}
// void GRDECLGrid::femTetCellWriter(const std::string &outputPath) const
// {
//     std::ofstream outputFile(outputPath);

//     if (!outputFile)
//     {
//         std::cerr << " Can not open: " << outputPath << std::endl;
//         std::abort();
//     }

//     outputFile << "NODES\n";
//     /**********Write Points Data******************/
//     for (size_t i = 0; i < allTetCoord_.size(); i++)
//     {
//         outputFile << allTetCoord_[i].x << ", ";
//         outputFile << allTetCoord_[i].y << ", ";
//         outputFile << allTetCoord_[i].z << "\n";
//     }
//     outputFile << "\n";
//     /**********Write CEll Data************/
//     // CELL header
//     // size_t cellNodesCounter = 0;
//     // for (size_t i = 0; i < decomposedCellNodes_.size(); i++)
//     // {
//     //     cellNodesCounter += decomposedCellNodes_[i].size();
//     // }
//     // int sizeTypeMap[8] = {-1, -1, -1, 10, 14, -1, -1, 12};
//     // Writing element data;
//     outputFile << "ELEMENTS\n";
//     for (size_t i = 0; i < allTetCellNodes_.size(); i++)
//     {
//         outputFile << i << ", ";
//         for (size_t j = 0; j < allTetCellNodes_[i].size(); j++)
//         {
//             outputFile << allTetCellNodes_[i][j];
//             if (j != allTetCellNodes_[i].size() - 1)
//             {
//                 outputFile << ", ";
//             }
//         }
//         outputFile << "\n";
//     }
//     outputFile << "\n";
// };
// void GRDECLGrid::femHexCellWriter(const std::string &outputPath) const
// {
//     std::ofstream outputFile(outputPath);

//     if (!outputFile)
//     {
//         std::cerr << " Can not open: " << outputPath << std::endl;
//         std::abort();
//     }

//     outputFile << "NODES\n";
//     /**********Write Points Data******************/
//     for (size_t i = 0; i < allHexCoord_.size(); i++)
//     {
//         outputFile << allHexCoord_[i].x << ", ";
//         outputFile << allHexCoord_[i].y << ", ";
//         outputFile << allHexCoord_[i].z << "\n";
//     }
//     outputFile << "\n";
//     /**********Write CEll Data************/
//     // CELL header
//     // size_t cellNodesCounter = 0;
//     // for (size_t i = 0; i < decomposedCellNodes_.size(); i++)
//     // {
//     //     cellNodesCounter += decomposedCellNodes_[i].size();
//     // }
//     // int sizeTypeMap[8] = {-1, -1, -1, 10, 14, -1, -1, 12};
//     // Writing element data;
//     outputFile << "ELEMENTS\n";
//     for (size_t i = 0; i < allHexCellNodes_.size(); i++)
//     {
//         outputFile << i << ", ";
//         for (size_t j = 0; j < allHexCellNodes_[i].size(); j++)
//         {
//             outputFile << allHexCellNodes_[i][j];
//             if (j != allHexCellNodes_[i].size() - 1)
//             {
//                 outputFile << ", ";
//             }
//         }
//         outputFile << "\n";
//     }
//     outputFile << "\n";
// };
void GRDECLGrid::writeRequiredData(const std::string& outputPath)
{
    DoubleList globalYoungMod(rawGRDECL_.dimI * rawGRDECL_.dimJ * rawGRDECL_.dimK, 0.1);
    DoubleList globalPoissonR(rawGRDECL_.dimI * rawGRDECL_.dimJ * rawGRDECL_.dimK, 0.2);
    if(localCellGlobalIndex_.size()!= YoungsMod_.size())
    {
        std::cerr << "Error: localCellGlobalIndex_ size not match with YoungsMod_!" << std::endl;
        std::abort();
    }
    for(size_t i = 0; i < localCellGlobalIndex_.size(); i++)
    {
        globalYoungMod[localCellGlobalIndex_[i]] = YoungsMod_[i];
        globalPoissonR[localCellGlobalIndex_[i]] = Poisson_[i];
    }
    // write to file
    std::ofstream outputFile(outputPath);
    if (!outputFile)
    {
        std::cerr << " Can not open: " << outputPath << std::endl;
        std::abort();
    }
    // write YOUNGMOD
    outputFile << std::setprecision(11);
    outputFile << "YOUNGMOD\n";
    for (size_t i = 0; i < globalYoungMod.size(); i++)
    {
        outputFile <<"1*"<< globalYoungMod[i] << " ";
    }
    outputFile << "\n";
    outputFile << "/";
    outputFile << "\n";
    // write POISSONR
    outputFile << "POISSONR\n";
    for (size_t i = 0; i < globalPoissonR.size(); i++)
    {
        outputFile <<"1*"<< globalPoissonR[i] << " ";
    }
    outputFile << "\n";
    outputFile << "/";
    outputFile << "\n";

}

void GRDECLGrid::vtkBoundaryWriter(const std::string &outputPath) const
{
    std::ofstream outputFile(outputPath);

    if (!outputFile)
    {
        std::cerr << " Can not open: " << outputPath << std::endl;
        std::abort();
    }

    /**********Prepare Data***********/
    Map<int, int> addedPointsMap; // Also: global to loacl Map
    IntList localPointToGlobalMap;
    for(IntSet::iterator it = KMinFace_.begin(); it != KMinFace_.end(); it++)
    {
        const int faceIdx = *it;
        const Face &faceI = faces_[faceIdx];
        for (size_t ni = 0; ni < faceI.size(); ni++)
        {
            const int nodeIdx = faceI[ni];
            if (addedPointsMap.find(nodeIdx) == addedPointsMap.end())
            {
                addedPointsMap.emplace(nodeIdx, int(localPointToGlobalMap.size()));
                localPointToGlobalMap.push_back(nodeIdx);
            }
        }
    }

    /**********Write header***********/
    outputFile << "# vtk DataFile Version 2.0\n";
    outputFile << "GRDECL Mesh\n";
    outputFile << "ASCII\n";
    outputFile << "DATASET UNSTRUCTURED_GRID\n";
    outputFile << "POINTS " << localPointToGlobalMap.size() << " double"
               << "\n";
    /**********Write Points Data******************/
    outputFile << std::setprecision(11);
    for (size_t i = 0; i < localPointToGlobalMap.size(); i++)
    {
        const int gIdx = localPointToGlobalMap[i];
        outputFile << nodeCoord_[gIdx].x << "  ";
        outputFile << nodeCoord_[gIdx].y << "  ";
        outputFile << nodeCoord_[gIdx].z << "\n";
    }
    outputFile << "\n";
    int counter = 0;
    for (IntSet::iterator it = KMinFace_.begin(); it != KMinFace_.end(); it++)
    {
        counter += int(faces_[*it].size());
    }

    outputFile << "CELLS " << KMinFace_.size() << " " << KMinFace_.size() + counter << "\n";
    for (IntSet::iterator it = KMinFace_.begin(); it != KMinFace_.end(); it++)
    {
        const int faceIdx = *it;
        outputFile << faces_[faceIdx].size() << " ";
        for (size_t j = 0; j < faces_[faceIdx].size(); j++)
        {
            outputFile << addedPointsMap[faces_[faceIdx][j]] << " ";
        }
        outputFile << "\n";
    }
    outputFile << "\n";

    outputFile << "CELL_TYPES " << KMinFace_.size() << '\n';
    for (size_t i = 0; i < KMinFace_.size(); i++)
    {
        outputFile << "9\n";
    }
    outputFile << '\n';    
}

void GRDECLGrid::vtkKMaxWriter(const std::string &outputPath) const
{
    std::ofstream outputFile(outputPath);
    if (!outputFile)
    {
        std::cerr << " Can not open: " << outputPath << std::endl;
        std::abort();
    }
    /**********Prepare Data***********/
    Map<int, int> addedPointsMap; // Also: global to loacl Map
    IntList localPointToGlobalMap;
    for (size_t i = 0; i<KMaxFaceNodes_.size(); i++)
    {
        const Face &faceI = KMaxFaceNodes_[i];
        for (size_t ni = 0; ni < faceI.size(); ni++)
        {
            const int nodeIdx = faceI[ni];
            if (addedPointsMap.find(nodeIdx) == addedPointsMap.end())
            {
                addedPointsMap.emplace(nodeIdx, int(localPointToGlobalMap.size()));
                localPointToGlobalMap.push_back(nodeIdx);
            }
        }
    }

    /**********Write header***********/
    outputFile << "# vtk DataFile Version 2.0\n";
    outputFile << "GRDECL Mesh\n";
    outputFile << "ASCII\n";
    outputFile << "DATASET UNSTRUCTURED_GRID\n";
    outputFile << "POINTS " << localPointToGlobalMap.size() << " double"
               << "\n";
    /**********Write Points Data******************/
    outputFile << std::setprecision(11);
    for (size_t i = 0; i < localPointToGlobalMap.size(); i++)
    {
        const int gIdx = localPointToGlobalMap[i];
        outputFile << allHexCoord_[gIdx].x << "  ";
        outputFile << allHexCoord_[gIdx].y << "  ";
        outputFile << allHexCoord_[gIdx].z << "\n";
    }
    outputFile << "\n";
    int counter = 0;
    for (size_t i = 0; i<KMaxFaceNodes_.size(); i++)
    {
        counter += int(KMaxFaceNodes_[i].size());
    }

    outputFile << "CELLS " << KMaxFaceNodes_.size() << " " << KMaxFaceNodes_.size() + counter << "\n";
    for (size_t i = 0; i<KMaxFaceNodes_.size(); i++)
    {
        const Face &faceI = KMaxFaceNodes_[i];
        outputFile << faceI.size() << " ";
        for (size_t j = 0; j < faceI.size(); j++)
        {
            outputFile << addedPointsMap[faceI[j]] << " ";
        }
        outputFile << "\n";
    }
    outputFile << "\n";

    outputFile << "CELL_TYPES " << KMaxFaceNodes_.size() << '\n';
    for (size_t i = 0; i < KMaxFaceNodes_.size(); i++)
    {
        outputFile << "9\n";
    }
    outputFile << '\n';        
}
void GRDECLGrid::vtkExtendedCell(const std::string &outputPath)
{
    // std::ofstream outputFile(outputPath);

    // if (!outputFile)
    // {
    //     std::cerr << " Can not open: " << outputPath << std::endl;
    //     std::abort();
    // }
    // /**********Write header***********/
    // outputFile << "# vtk DataFile Version 2.0\n";
    // outputFile << "GRDECL Mesh\n";
    // outputFile << "ASCII\n";
    // outputFile << "DATASET UNSTRUCTURED_GRID\n";
    // outputFile << "POINTS " << allHexCoord_.size() << " double"
    //            << "\n";
    // /**********Write Points Data******************/
    // outputFile << std::setprecision(11);
    // for (size_t i = 0; i < allHexCoord_.size(); i++)
    // {
    //     outputFile << allHexCoord_[i].x << "  ";
    //     outputFile << allHexCoord_[i].y << "  ";
    //     outputFile << allHexCoord_[i].z << "\n";
    // }
    // outputFile << "\n";
    // /**********Write CEll Data************/
    // // CELL header
    // int totalCellNum = int(allHexCellNodes_.size() + overBurdenCellNodes_.size() + underBurdenCellNodes_.size() + IMaxFaceCellNodes_.size() + IMinFaceCellNodes_.size() + JMaxFaceCellNodes_.size()+JMinFaceCellNodes_.size());
    // outputFile << "CELLS " << totalCellNum << " " <<totalCellNum * 9 << "\n";
    // for (size_t i = 0; i < allHexCellNodes_.size(); i++)
    // {
    //     outputFile << "8 ";
    //     for (size_t j = 0; j < allHexCellNodes_[i].size(); j++)
    //     {
    //         outputFile << allHexCellNodes_[i][j] << " ";
    //     }
    //     outputFile << "\n";
    // }
    // for (size_t i = 0; i < overBurdenCellNodes_.size(); i++)
    // {
    //     outputFile << "8 ";
    //     for (size_t j = 0; j < overBurdenCellNodes_[i].size(); j++)
    //     {
    //         outputFile << overBurdenCellNodes_[i][j] << " ";
    //     }
    //     outputFile << "\n";
    // }
    // for (size_t i = 0; i < underBurdenCellNodes_.size(); i++)
    // {
    //     outputFile << "8 ";
    //     for (size_t j = 0; j < underBurdenCellNodes_[i].size(); j++)
    //     {
    //         outputFile << underBurdenCellNodes_[i][j] << " ";
    //     }
    //     outputFile << "\n";
    // }
    // outputFile << "\n";
    // for (size_t i = 0; i < IMaxFaceCellNodes_.size(); i++)
    // {
    //     outputFile << "8 ";
    //     for (size_t j = 0; j < IMaxFaceCellNodes_[i].size(); j++)
    //     {
    //         outputFile << IMaxFaceCellNodes_[i][j] << " ";
    //     }
    //     outputFile << "\n";
    // }
    // for (size_t i = 0; i < IMinFaceCellNodes_.size(); i++)
    // {
    //     outputFile << "8 ";
    //     for (size_t j = 0; j < IMinFaceCellNodes_[i].size(); j++)
    //     {
    //         outputFile << IMinFaceCellNodes_[i][j] << " ";
    //     }
    //     outputFile << "\n";
    // }
    // for(size_t i = 0; i < JMaxFaceCellNodes_.size(); i++)
    // {
    //     outputFile << "8 ";
    //     for (size_t j = 0; j < JMaxFaceCellNodes_[i].size(); j++)
    //     {
    //         outputFile << JMaxFaceCellNodes_[i][j] << " ";
    //     }
    //     outputFile << "\n";
    // }
    // for(size_t i = 0; i < JMinFaceCellNodes_.size(); i++)
    // {
    //     outputFile << "8 ";
    //     for (size_t j = 0; j < JMinFaceCellNodes_[i].size(); j++)
    //     {
    //         outputFile << JMinFaceCellNodes_[i][j] << " ";
    //     }
    //     outputFile << "\n";
    // }

    // // Writing cell types;
    // outputFile << "CELL_TYPES " << totalCellNum << '\n';
    // for (size_t i = 0; i < totalCellNum; i++)
    // {
    //     outputFile << "12\n";
    // }

    // /**********Write POINT Data************/
    // //prepare data 
    // IntSet faultPoints; // Also: global to loacl Map
    // for (size_t fi = 0; fi < faultFaces_.size(); fi++)
    // {
    //     const int faceIdx = faultFaces_[fi];
    //     const Face &faceI = faces_[faceIdx];
    //     for (size_t ni = 0; ni < faceI.size(); ni++)
    //     {
    //         const int nodeIdx = faceI[ni];
    //         if (faultPoints.find(nodeIdx) == faultPoints.end())
    //         {
    //             faultPoints.insert(nodeIdx);
    //         }
    //     }
    // }
    // outputFile << "POINT_DATA " << allHexCoord_.size() << '\n';
    // outputFile << "SCALARS FaultTag double 1\n";
    // outputFile << "LOOKUP_TABLE FaultTag\n";
    // for (size_t i = 0; i < allHexCoord_.size(); i++)
    // {
    //     if (faultPoints.find(int(i)) != faultPoints.end())
    //     {
    //         outputFile << "1\n";
    //     }
    //     else
    //     {
    //         outputFile << "0\n";
    //     }
    // }

    // outputFile << '\n';

    // // free data for memory test
    // allHexCoord_.resize(0);
    // allHexCellNodes_.resize(0);
}

void GRDECLGrid::vtkExtendedCell(
    const std::string &outputPath,
    const PointList &nodeCoord,
    const IntListList &cellNodes,
    const IntList &pointData,
    const IntList &cellData)
{
    std::ofstream outputFile(outputPath);
    if (!outputFile)
    {
        std::cerr << " Can not open: " << outputPath << std::endl;
        std::abort();
    }
    if (nodeCoord.size() != pointData.size())
    {
        std::cerr << "Error: pointData size not match with nodeCoord size!" << std::endl;
        std::abort();
    }
    if (cellNodes.size() != cellData.size())
    {
        std::cerr << "Error: cellData size not match with cellNodes size!" << std::endl;
        std::abort();
    }

    /**********Write header***********/
    outputFile << "# vtk DataFile Version 2.0\n";
    outputFile << "GRDECL Mesh\n";
    outputFile << "ASCII\n";
    outputFile << "DATASET UNSTRUCTURED_GRID\n";
    outputFile << "POINTS " << nodeCoord.size() << " double"
               << "\n";
    /**********Write Points Data******************/
    outputFile << std::setprecision(11);
    for (size_t i = 0; i < nodeCoord.size(); i++)
    {
        outputFile << nodeCoord[i].x << "  ";
        outputFile << nodeCoord[i].y << "  ";
        outputFile << nodeCoord[i].z << "\n";
    }
    outputFile << "\n";
    /**********Write CEll Data************/
    // CELL header
    // loop cell for find sum
    size_t cellNodesCounter = 0;
    for (size_t i = 0; i < cellNodes.size(); i++)
    {
        cellNodesCounter += cellNodes[i].size();
    }
    int sizeTypeMap[8] = {-1, -1, -1, 10, 14, 13, -1, 12};
    // Writing cell data;
    outputFile << "CELLS " << cellNodes.size() << " " << cellNodes.size() + cellNodesCounter << "\n";
    for (size_t i = 0; i < cellNodes.size(); i++)
    {
        outputFile << cellNodes[i].size() << " ";
        for (size_t j = 0; j < cellNodes[i].size(); j++)
        {
            outputFile << cellNodes[i][j] << " ";
        }
        outputFile << "\n";
    }
    outputFile << "\n";

    // Writing cell types;
    outputFile << "CELL_TYPES " << cellNodes.size() << '\n';
    for (size_t i = 0; i < cellNodes.size(); i++)
    {
        outputFile << sizeTypeMap[cellNodes[i].size() - 1] << '\n';
    }

    /**********Write POINT Data************/
    outputFile << "POINT_DATA " << pointData.size() << '\n';
    outputFile << "SCALARS PointTag int 1\n";
    outputFile << "LOOKUP_TABLE PointTag\n";
    for (size_t i = 0; i < pointData.size(); i++)
    {
        outputFile << pointData[i] << "\n";
    }
    outputFile << '\n';
    /**********Write POINT Data************/
    outputFile << "CELL_DATA " << cellNodes.size() << '\n';
    outputFile << "SCALARS CellTag int 1\n";
    outputFile << "LOOKUP_TABLE CellTag\n";
    for (size_t i = 0; i < cellData.size(); i++)
    {
        outputFile << cellData[i] << "\n";
    }
    outputFile << '\n';
}

void GRDECLGrid::genHybridExtendCell() {
    
    genExtendCell(decomposedCoord_,decomposedCellNodes_,decomposedPointData_,decomposedCellData_);
    decomposeHybridExtendCell();
    //
    decomposedYoungsMod_.resize(decomposedCellNodes_.size(),0.0);
    decomposedPoisson_.resize(decomposedCellNodes_.size(),0.0);
    decomposedPoro_.resize(decomposedCellNodes_.size(),0.0);
    decomposedRockDens_.resize(decomposedCellNodes_.size(),0.0);
    //
    decomposedCellVolume_.resize(decomposedCellNodes_.size(),0.0);
    decomposedCellCenter_.resize(decomposedCellNodes_.size());
    calcExtendCellCenterAndVolume(decomposedCoord_,decomposedCellNodes_,decomposedCellCenter_,decomposedCellVolume_);
}
void GRDECLGrid::genHexExtendCell() {
    genExtendCell(allHexCoord_,allHexCellNodes_,hexPointData_,hexCellData_);
    //
    hexYoungsMod_.resize(allHexCellNodes_.size(),0.0);
    hexPoisson_.resize(allHexCellNodes_.size(),0.0);
    hexPoro_.resize(allHexCellNodes_.size(),0.0);
    hexRockDens_.resize(allHexCellNodes_.size(),0.0);
    // 
    hexCellVolume_.resize(allHexCellNodes_.size(),0.0);
    hexCellCenter_.resize(allHexCellNodes_.size());
    calcExtendCellCenterAndVolume(allHexCoord_,allHexCellNodes_,hexCellCenter_,hexCellVolume_);
}
void GRDECLGrid::genTetExtendCell() {
    genExtendCell(allTetCoord_,allTetCellNodes_,tetPointData_,tetCellData_);
    decomposeAllTetExtendCell();
    //
    tetYoungsMod_.resize(allTetCellNodes_.size(),0.0);
    tetPoisson_.resize(allTetCellNodes_.size(),0.0);
    tetPoro_.resize(allTetCellNodes_.size(),0.0);
    tetRockDens_.resize(allTetCellNodes_.size(),0.0);
    //
    tetCellVolume_.resize(allTetCellNodes_.size(),0.0);
    tetCellCenter_.resize(allTetCellNodes_.size());
    calcExtendCellCenterAndVolume(allTetCoord_,allTetCellNodes_,tetCellCenter_,tetCellVolume_);
}

void GRDECLGrid::vtkExtendedCellData(
    const std::string &outputPath,
    const PointList &nodeCoord,
    const IntListList &cellNodes,
    const IntList &pointData,
    const IntList &cellData,
    const DoubleList &youngsMod,
    const DoubleList &poisson,
    const DoubleList &poro,
    const DoubleList &rockDens)
{
    std::ofstream outputFile(outputPath);

    if (!outputFile)
    {
        std::cerr << " Can not open: " << outputPath << std::endl;
        std::abort();
    }
    if(nodeCoord.size() != pointData.size())
    {
        std::cerr << "Error: pointData size not match with nodeCoord size!" << std::endl;
        std::abort();
    }
    if (cellNodes.size() != cellData.size() ||
        cellNodes.size() != youngsMod.size() ||
        cellNodes.size() != poisson.size() ||
        cellNodes.size() != poro.size() ||
        cellNodes.size() != rockDens.size())
    {
        std::cerr << "Error: cellData size not match with cellNodes size!" << std::endl;
        std::abort();
    }

        /**********Write header***********/
    outputFile << "# vtk DataFile Version 2.0\n";
    outputFile << "GRDECL Mesh\n";
    outputFile << "ASCII\n";
    outputFile << "DATASET UNSTRUCTURED_GRID\n";
    outputFile << "POINTS " << nodeCoord.size() << " double"
               << "\n";
    /**********Write Points Data******************/
    outputFile << std::setprecision(11);
    for (size_t i = 0; i < nodeCoord.size(); i++)
    {
        outputFile << nodeCoord[i].x << "  ";
        outputFile << nodeCoord[i].y << "  ";
        outputFile << nodeCoord[i].z << "\n";
    }
    outputFile << "\n";
    /**********Write CEll Data************/
    // CELL header
    // loop cell for find sum
    size_t cellNodesCounter = 0;
    for (size_t i = 0; i < cellNodes.size(); i++)
    {
        cellNodesCounter += cellNodes[i].size();
    }
    int sizeTypeMap[8] = {-1, -1, -1, 10, 14, 13, -1, 12};
    // Writing cell data;
    outputFile << "CELLS " << cellNodes.size() << " " << cellNodes.size() + cellNodesCounter << "\n";
    for (size_t i = 0; i < cellNodes.size(); i++)
    {
        outputFile << cellNodes[i].size() << " ";
        for (size_t j = 0; j < cellNodes[i].size(); j++)
        {
            outputFile << cellNodes[i][j] << " ";
        }
        outputFile << "\n";
    }
    outputFile << "\n";

    // Writing cell types;
    outputFile << "CELL_TYPES " << cellNodes.size() << '\n';
    for (size_t i = 0; i < cellNodes.size(); i++)
    {
        outputFile << sizeTypeMap[cellNodes[i].size()-1] << '\n';
    }

    /**********Write POINT Data************/
    outputFile << "POINT_DATA " << pointData.size() << '\n';
    outputFile << "SCALARS PointTag int 1\n";
    outputFile << "LOOKUP_TABLE default\n";
    for (size_t i = 0; i < pointData.size(); i++)
    {
        outputFile << pointData[i] << "\n";
    }
    outputFile << '\n';
    /**********Write POINT Data************/
    outputFile << "CELL_DATA " << cellNodes.size() << '\n';
    outputFile << "SCALARS CellTag int 1\n";
    outputFile << "LOOKUP_TABLE default\n";
    for (size_t i = 0; i < cellData.size(); i++)
    {
        outputFile << cellData[i] << "\n";
    }
    outputFile << '\n';
    /**********Write pros Data************/
    outputFile << "SCALARS YoungModulos double 1\n";
    outputFile << "LOOKUP_TABLE default\n";
    for (size_t i = 0; i < youngsMod.size(); i++)
    {
        outputFile << youngsMod[i] << "\n";
    }
    outputFile << '\n';
    outputFile << "SCALARS PoissonRatio double 1\n";
    outputFile << "LOOKUP_TABLE default\n";
    for (size_t i = 0; i < poisson.size(); i++)
    {
        outputFile << poisson[i] << "\n";
    }
    outputFile << '\n';
    outputFile << "SCALARS RockDensity double 1\n";
    outputFile << "LOOKUP_TABLE default\n";
    for (size_t i = 0; i < rockDens.size(); i++)
    {
        outputFile << rockDens[i] << "\n";
    }
    outputFile << '\n';
    outputFile << "SCALARS Porosity double 1\n";
    outputFile << "LOOKUP_TABLE default\n";
    for (size_t i = 0; i < poro.size(); i++)
    {
        outputFile << poro[i] << "\n";
    }
    outputFile << '\n';
}

void GRDECLGrid::scaleCellData(DoubleList& data, const double scale) const
{
    for(size_t i = 0; i < data.size(); i++)
    {
        data[i] *= scale;
    }
}

void GRDECLGrid::assignCellData(
    const IntList &cellTypes,
    DoubleList &data,
    const double overValue,
    const double underValue,
    const double sideValue,
    const double overSideValue,
    const double underSideValue)
{   
    if (cellTypes.size() != data.size())
    {
        std::cerr << "Error: cellTypes size not match with data size!" << std::endl;
        std::abort();
    }
    for (size_t i = 0; i < cellTypes.size(); i++)
    {
        switch (cellTypes[i])
        {
        case 2:
            data[i] = sideValue;
            break;
        case 3:
            data[i] = overValue;
            break;
        case 4:
            data[i] = overSideValue;
            break;
        case 5:
            data[i] = underValue;
            break;
        case 6:
            data[i] = underSideValue;
            break;
        default:
            break;
        }
    }
}

void GRDECLGrid::calcExtendCellCenterAndVolume(
    const PointList &nodeCoord,
    const IntListList &cellNodes,
    PointList &cellCenter,
    DoubleList &cellVolume) const
{
    cellCenter.resize(cellNodes.size());
    cellVolume.resize(cellNodes.size());
    for (size_t i = 0; i < cellNodes.size(); i++)
    {
        Point center(0.0, 0.0, 0.0);
        cellCenter[i] = center;
        for (size_t j = 0; j < cellNodes[i].size(); j++)
        {
            center += nodeCoord[cellNodes[i][j]];
        }
        center /= double(cellNodes[i].size());
        cellCenter[i] = center;
        cellVolume[i] = calcCellVolumeByNodes(nodeCoord,cellCenter[i],cellNodes[i]);
    }
}

// 计算守恒值
double GRDECLGrid::calcConservation(const DoubleList &value, const DoubleList &coeff, const size_t dataSize)
{
    double conservation = 0;
    for (size_t i = 0; i < dataSize; i++)
    {
        conservation += value[i] * coeff[i];
    }
    return conservation;
}
// 缩放数据
void GRDECLGrid::scaleData(DoubleList &data, const size_t dataSize, const double scale)
{
    for (size_t i = 0; i < dataSize; i++)
    {
        data[i] *= scale;
    }
}
// 提取界面
void GRDECLGrid::extractLayer(IntList &pointData, const IntList &layers)
{
    Map<int, int> layerMap;
    for (size_t i = 0; i < layers.size(); i++)
    {
        layerMap.emplace(layers[i], int(i));
    }
    for (size_t i = 0; i < pointData.size(); i++)
    {
        if (pointData[i] > 9)
        {
            if (layerMap.find(pointData[i]) != layerMap.end())
            {
                pointData[i] = layerMap[pointData[i]] + 10;
            }
            else
            {
                pointData[i] = 0;
            }
        }
    }
}

void GRDECLGrid::freeGrdeclData()
{
    //localCellGlobalIndex_.resize(0);
    //globalCellLocalIndex_.resize(0);
    cornerMap_.resize(0);
    posPointIndex_.resize(0);
    intersectIndices_.resize(0);
    //cellVolumes_.resize(0);
    //faceCenters_.resize(0);
    //faceAreaNormals_.resize(0);
    //cellNodes_.resize(0);
    //cellFaces_.resize(0);
    //faultFaces_.resize(0);
    //faceTags_.resize(0);
    //faceNeighbors_.resize(0);
    //nodeCoord_.resize(0);
    //faces_.resize(0);
}
void GRDECLGrid::freeHexData()
{
    allHexCoord_.resize(0);
    allHexCellNodes_.resize(0);
    hexPointData_.resize(0);
    hexCellData_.resize(0);
    hexYoungsMod_.resize(0);
    hexPoisson_.resize(0);
    hexPoro_.resize(0);
    hexRockDens_.resize(0);
    hexCellVolume_.resize(0);
    hexCellCenter_.resize(0);
}

void GRDECLGrid::freeTetData()
{
    allTetCoord_.resize(0);
    allTetCellNodes_.resize(0);
    tetPointData_.resize(0);
    tetCellData_.resize(0);
    tetYoungsMod_.resize(0);
    tetPoisson_.resize(0);
    tetPoro_.resize(0);
    tetRockDens_.resize(0);
    tetCellVolume_.resize(0);
    tetCellCenter_.resize(0);
}

void GRDECLGrid::freeDecomposedData()
{
    decomposedCoord_.resize(0);
    decomposedCellNodes_.resize(0);
    decomposedPointData_.resize(0);
    decomposedCellData_.resize(0);
    decomposedYoungsMod_.resize(0);
    decomposedPoisson_.resize(0);
    decomposedPoro_.resize(0);
    decomposedRockDens_.resize(0);
    decomposedCellVolume_.resize(0);
    decomposedCellCenter_.resize(0);
}